DOI

Project overview

The goal of the Pyromania project is to test how terrestrial subsides (plant biomass loading or “browning”) and burning influence aquatic productivity, water quality/chemistry, and trophic transfer. We used a manipulative experiment to assess a range of plant material quantities (0-400g per tank) and fire treatment (burned vs unburned material) and the non-linearity of these effects on aquatic systems through 4 time-point samplings. We used 400L aquatic mesocosms and ran the experiment for ~90d in 2021-2022.

Figure. Images showing (A) controlled burning of plant material, (B) dried unburned plants and (C) packed plant material in mesh bags, (D) visual differences in water quality between tanks receiving high (left) and low (right) amounts of plant detritus, and (E) the experimental mesocosm array at the University of San Diego Biological Field Station.
Figure. Images showing (A) controlled burning of plant material, (B) dried unburned plants and (C) packed plant material in mesh bags, (D) visual differences in water quality between tanks receiving high (left) and low (right) amounts of plant detritus, and (E) the experimental mesocosm array at the University of San Diego Biological Field Station.

DATA SETS
This data set is among 3 to be generated for the project and focuses on:

  1. dissolved organic carbon (DOC) concentrations
  2. total dissolved nitrogen (TDN) concentrations
  3. net primary productivity and respiration (NPP, Resp) using whole mesocosm dissolved oxygen (DO) measurements
  4. burning effects on donor plant material using elemental analysis
  5. trophic transfer using 15N-labeled sage plants to trace N from the plants into plankton.
  6. greenhouse gas concentrations (carbon dioxide (CO2), methane (CH4)) from mesocosms

TIME POINTS

  • Time-0 (T0) = before the addition of plant materials, plankton were stocked at this stage
  • Time-1 (T1) = 10 days after the addition of plant materials
  • Time-2 (T2) = 31 days after …
  • Time-3 (T3) = 59 days after …
  • Time-4 (T4) = 89 days after …

General notes on GAM analyses
We fit the generalized additive models (GAMs) via restricted maximum likelihood (REML) to give stable results with the smoothing parameter (sp) to determine the non-linear relationship between response variables and plant-biomass loading (x-axis). We use automatic smoothing with k value generated automatically from the models, which will set the line ‘wiggliness’. Too low and the relationship becomes linear; too high, and the wiggliness goes haywire.

When using the non-linear smoothing, this is the s(x). When the variable is inside the smooth function, this accounts for the nonlinear shape. We do not use additive non-linear smoothing, which is when two smoothers together, as s(x1) + s(x2), instead we use factor-smooth interaction (detailed below). In addition, we use Treatment (and occasionally plankton size fractions, or Type) as predictors outside of the smooth terms s(x1); this allows for linearity. Continuous variables are rarely linear in GAMs, however, setting categorical variables as linear predictors is more common.

Factor-smooth interactions are written as s(x1 by = fac). This sets different smoothers for different variables of “fac”. Usually, the different smoothers are combined with a varying intercept in case the different categories are different in means and slopes, this would be by adding the fac + s(x1 by = fac), where the +fac allows for a different slope. Similarly, in the absence of by = fac, the smoother is considered a global smoother s(x1), fitting a single line to all the data. If a global smoother is combined with a factor term, then this is akin to varying the intercept but keeping the same slope: fac + s(x1).

The EDF - effective degrees of freedom equate with wiggliness, where edf =1 is a straight line, and higher edfs as more wiggly. GAM smoother significance is described as not being able to draw a horizontal line through the data. Finally, it is also advised to check model concurvity, which is the collinearity with models from 0-1.

Plant material

Plant material for the starting material (sage or willow, stems or sticks). This is useful in determining how fire impacted the nutrients in the plant material.

  • First, run some stats to see what is happening and where differences lie.
  • We will then make a summary boxplot figure
plant.nut<-read.csv("data/Pyro_plant material_elemental.csv")

plant.N<-aggregate(N~plant+treatment, plant.nut, FUN=mean)

fac<-c("Sample.Name", "type", "plant", "treatment", "rep") # make factors
plant.nut[fac]<-lapply(plant.nut[fac],factor)

######### all plant N test
plant.N<-lm(N~treatment+type, data=plant.nut)
Anova(plant.N, type=3) # 2 way interaction and main effects
## Anova Table (Type III tests)
## 
## Response: N
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 21.0899  1 611.5406 < 2.2e-16 ***
## treatment    0.0410  1   1.1896    0.2884    
## type         3.2033  1  92.8850 5.863e-09 ***
## Residuals    0.6897 20                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# use the above to address differences in %N for sage-derived 15N x-axis = plant N added to tanks

# look at mean, SE, n for the %N
plant.mean.N<-aggregate(N~treatment, plant.nut, FUN=mean) 
plant.SE.N<-aggregate(N~treatment, plant.nut, FUN=std.error)
plant.n<-aggregate(N~treatment, plant.nut, FUN=length)

####### separate dfs by species and test elemental analysis

sage.nut<-plant.nut[(plant.nut$plant=="sage"),]
will.nut<-plant.nut[(plant.nut$plant=="willow"),]

### Sage ### 
######### %N
plant.N.sag<-lm(N~treatment*type, data=sage.nut)
Anova(plant.N.sag, type=3) # 2 way interaction and main effects
## Anova Table (Type III tests)
## 
## Response: N
##                Sum Sq Df   F value    Pr(>F)    
## (Intercept)    7.0227  1 2171.5213 4.971e-11 ***
## treatment      0.0280  1    8.6632  0.018615 *  
## type           0.2604  1   80.5246 1.894e-05 ***
## treatment:type 0.0705  1   21.8099  0.001602 ** 
## Residuals      0.0259  8                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.N.sag, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean     SE df lower.CL upper.CL .group
##  burned     1.530 0.0328  8    1.454     1.61  a    
##  unburned   1.667 0.0328  8    1.591     1.74   b   
## 
## type = stem:
##  treatment emmean     SE df lower.CL upper.CL .group
##  unburned   0.943 0.0328  8    0.868     1.02  a    
##  burned     1.113 0.0328  8    1.038     1.19   b   
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

plant.part.N<-aggregate(N~treatment+type, sage.nut, FUN=mean)
#burning reduced Sage leaf N = (1.53-1.67)/1.67  = -8%
#burning increased Sage stem N= (1.11-0.94)/0.94  = 18%

########## %P 
plant.P.sag<-lm(P~treatment*type, data=sage.nut)
Anova(plant.P.sag, type=3) # just type
## Anova Table (Type III tests)
## 
## Response: P
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.314928  1 804.0715 2.586e-09 ***
## treatment      0.000771  1   1.9677  0.198292    
## type           0.007561  1  19.3060  0.002306 ** 
## treatment:type 0.000056  1   0.1438  0.714371    
## Residuals      0.003133  8                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


########## %K
plant.K.sag<-lm(K~treatment*type, data=sage.nut)
Anova(plant.K.sag, type=3) # type and treatment
## Anova Table (Type III tests)
## 
## Response: K
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)    12.5256  1 372.2328 5.404e-08 ***
## treatment       0.2563  1   7.6157   0.02469 *  
## type            0.1803  1   5.3571   0.04934 *  
## treatment:type  0.0736  1   2.1882   0.17733    
## Residuals       0.2692  8                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.K.sag, ~treatment)
multcomp::cld(posthoc, Letters=letters)
##  treatment emmean     SE df lower.CL upper.CL .group
##  unburned    1.61 0.0749  8     1.44     1.79  a    
##  burned      1.87 0.0749  8     1.70     2.04   b   
## 
## Results are averaged over the levels of: type 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

plant.part.K<-aggregate(K~treatment, sage.nut, FUN=mean)
#burning increased Sage leaf K = (2.04-1.63)/1.63  = 25%

########## %S
plant.S.sag<-lm(S~treatment*type, data=sage.nut)
Anova(plant.S.sag, type=3) # type effect
## Anova Table (Type III tests)
## 
## Response: S
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.67783  1 859.3665 1.986e-09 ***
## treatment      0.00010  1   0.1321    0.7257    
## type           0.11152  1 141.3891 2.299e-06 ***
## treatment:type 0.00039  1   0.4885    0.5044    
## Residuals      0.00631  8                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


############## Zn ppm
plant.ZN.sag<-lm(ZN~treatment*type, data=sage.nut)
Anova(plant.ZN.sag, type=3) #no effect
## Anova Table (Type III tests)
## 
## Response: ZN
##                 Sum Sq Df F value   Pr(>F)   
## (Intercept)    19959.4  1 15.1723 0.004576 **
## treatment       3073.6  1  2.3364 0.164903   
## type            3280.7  1  2.4938 0.152947   
## treatment:type  1236.3  1  0.9398 0.360731   
## Residuals      10524.1  8                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.ZN.sag, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    36.3 20.9  8    -12.0     84.6  a    
##  burned      81.6 20.9  8     33.3    129.9  a    
## 
## type = stem:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    30.1 20.9  8    -18.2     78.4  a    
##  burned      34.8 20.9  8    -13.5     83.1  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.


###############################################
###############################################
### Willow ### 
######### %N
plant.N.wil<-lm(N~treatment*type, data=will.nut)
Anova(plant.N.wil, type=3) #  main effects
## Anova Table (Type III tests)
## 
## Response: N
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept)    8.7381  1 644.0658 3.766e-08 ***
## treatment      0.0817  1   6.0194   0.04388 *  
## type           1.4538  1 107.1530 1.703e-05 ***
## treatment:type 0.0436  1   3.2119   0.11620    
## Residuals      0.0950  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plant.part.Nwill<-aggregate(N~treatment+type, will.nut, FUN=mean)
#burning increased Sage leaf K = (1.27-1.05)/1.05  = 21%


########## %P 
plant.P.wil<-lm(P~treatment*type, data=will.nut)
Anova(plant.P.wil, type=3) # type and treatment
## Anova Table (Type III tests)
## 
## Response: P
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.156865  1 447.9429 1.323e-07 ***
## treatment      0.006403  1  18.2834  0.003674 ** 
## type           0.007239  1  20.6703  0.002647 ** 
## treatment:type 0.000880  1   2.5131  0.156921    
## Residuals      0.002451  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.P.wil, ~treatment)
multcomp::cld(posthoc, Letters=letters)
##  treatment emmean      SE df lower.CL upper.CL .group
##  unburned   0.143 0.00764  7    0.125    0.161  a    
##  burned     0.190 0.00854  7    0.170    0.210   b   
## 
## Results are averaged over the levels of: type 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

plant.part.Pwill<-aggregate(P~treatment, will.nut, FUN=mean)
#burning increased Sage leaf P = (0.20-0.14)/0.14  = 43%


########## %K
plant.K.wil<-lm(K~treatment*type, data=will.nut)
Anova(plant.K.wil, type=3) # type and treatment
## Anova Table (Type III tests)
## 
## Response: K
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept)    4.6128  1 548.0120  6.59e-08 ***
## treatment      0.0676  1   8.0344 0.0252428 *  
## type           0.2640  1  31.3583 0.0008161 ***
## treatment:type 0.0015  1   0.1725 0.6903484    
## Residuals      0.0589  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.K.wil, ~treatment)
multcomp::cld(posthoc, Letters=letters)
##  treatment emmean     SE df lower.CL upper.CL .group
##  unburned   0.817 0.0375  7    0.728    0.905  a    
##  burned     1.006 0.0419  7    0.906    1.105   b   
## 
## Results are averaged over the levels of: type 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

plant.part.Kwill<-aggregate(K~treatment, will.nut, FUN=mean)
#burning increased Sage leaf K = (1.05-0.82)/0.82  = 28%


########## %S
plant.S.wil<-lm(S~treatment*type, data=will.nut)
Anova(plant.S.wil, type=3) # main and interactions
## Anova Table (Type III tests)
## 
## Response: S
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)    0.197633  1 1748.117 1.168e-09 ***
## treatment      0.002521  1   22.303  0.002151 ** 
## type           0.042375  1  374.819 2.446e-07 ***
## treatment:type 0.001355  1   11.985  0.010519 *  
## Residuals      0.000791  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.S.wil, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean      SE df lower.CL upper.CL .group
##  unburned  0.2157 0.00614  7   0.2012   0.2302  a    
##  burned    0.2567 0.00614  7   0.2422   0.2712   b   
## 
## type = stem:
##  treatment emmean      SE df lower.CL upper.CL .group
##  burned    0.0688 0.00752  7   0.0510   0.0865  a    
##  unburned  0.0728 0.00614  7   0.0583   0.0873  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

plant.part.Swill<-aggregate(S~treatment+type, will.nut, FUN=mean)
#burning increased Sage leaf S = (0.26-0.22)/0.22  = 18%


############## Zn ppm
plant.ZN.wil<-lm(ZN~treatment*type, data=will.nut)
Anova(plant.ZN.wil, type=3) #no effect
## Anova Table (Type III tests)
## 
## Response: ZN
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept)    124440  1 196.4564 2.229e-06 ***
## treatment       22363  1  35.3043 0.0005748 ***
## type            21956  1  34.6631 0.0006071 ***
## treatment:type   5340  1   8.4306 0.0228697 *  
## Residuals        4434  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc<-emmeans(plant.ZN.wil, ~treatment| type)
multcomp::cld(posthoc, Letters=letters)
## type = leaf:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    81.6 14.5  7    47.21    115.9  a    
##  burned     203.7 14.5  7   169.31    238.0   b   
## 
## type = stem:
##  treatment emmean   SE df lower.CL upper.CL .group
##  unburned    35.8 14.5  7     1.44     70.2  a    
##  burned      68.4 17.8  7    26.32    110.5  a    
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

plant.part.Znwill<-aggregate(ZN~treatment+type, will.nut, FUN=mean)
#burning increased Sage leaf K = (204-82)/82  = 148%
plant.nut$int.fac<-factor(interaction(plant.nut$plant, plant.nut$type),
                         levels=c("sage.leaf", "sage.stem", "willow.leaf", "willow.stem"))
####### figures
N.plot<- ggplot(plant.nut, aes(x=int.fac, y=N, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%N")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

P.plot<- ggplot(plant.nut, aes(x=int.fac, y=P, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%P")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

K.plot<- ggplot(plant.nut, aes(x=int.fac, y=K, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%K")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))


S.plot<- ggplot(plant.nut, aes(x=int.fac, y=S, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75))+
  xlab("Plant:Tissue")+
  ylab("%S")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Zn.plot<- ggplot(plant.nut, aes(x=int.fac, y=ZN, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
                 position=position_dodge(0.75)) +
  xlab("Plant:Tissue")+
  ylab("Zn ppm")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
  
extract.legend.nut <- get_legend(
  # create some space to the left of the legend
  Zn.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))

nutrients<-plot_grid(
  N.plot+ theme(legend.position = "none"),
  P.plot+ theme(legend.position = "none"),
  K.plot+ theme(legend.position = "none"), 
  S.plot+ theme(legend.position = "none"), 
  Zn.plot+ theme(legend.position = "none"),
  extract.legend.nut, 
  rel_widths = c(8,8,8,8,8,3), ncol=6)

nutrients
**Figure** Elemental analysis of burned and unburned plant material (leaves and stem) from sage and willow prior to being added to experimental treatments.

Figure Elemental analysis of burned and unburned plant material (leaves and stem) from sage and willow prior to being added to experimental treatments.

ggsave("figures/leaf.nutrients.pdf", height=4, width=12)

DOC

Import the data for DOC and TDN and do a loop to clean up all files and make stacked data in single df. This will take the raw data files, align metadata, and filter to make a new df for models.

  • Analyze DOC at each time point.
  • Run model selection and produce plots for each individual timepoint
  • Graph the pooled timepoints as a 5 panel figure.
detach("package:dplyr", unload = TRUE)
library(dplyr)

## import treatment IDs
IDs<-read.csv("data/treatment.IDs.csv")

##### grab files in a list
Total.DOC.files <- list.files(path="data/DOC.TN", pattern = "csv$", full.names = T)

##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="data/DOC.TN", pattern = "csv$", full.names = F))
############ formatting all data in for loop
  for(i in 1:length(Total.DOC.files))
    {
  data<-read.csv(Total.DOC.files[i], sep=",")
  data<-data[,c(1:3)] # only keep these columns
  data$File<-Total.DOC.files[i]
  colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
  data$Tank<- IDs$Tank
  data$Tank<-as.numeric(as.character(data$Tank)) # make the column of samples all numeric
  data <- data[!is.na(as.numeric(as.character(data$Tank))),] # remove all rows that aren't numeric/tanks
  data$Treatment<-IDs$Treatment
  data$plant.mass..g<-IDs$plant.mass..g
  make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
  }
########## this is the end of the loop

#see all dfs you've made, the above will be df matching their file names
# ls() 

#Combine files from loop to single df
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)

DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names

# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27) 

#alternatively
# remove the 9 letters ('^.) at start 
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File) 

# if equals DOC_T0_11052021 then, T0, if not then T1
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
         ifelse (DOC.df$File=="DOC_T1.csv", "T1",
           ifelse (DOC.df$File=="DOC_T2.csv", "T2", 
                   ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))

#rearrange
DOC.df<- DOC.df %>% 
  select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L) 

DOC.df$Treatment<-as.factor(DOC.df$Treatment)

######## T0 model
m1.DOC.T0<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.DOC.T0<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.DOC.T0<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.DOC.AIC<-AIC(m1.DOC.T0, m2.DOC.T0, m3.DOC.T0)
# best is smoother solo

summary(m3.DOC.T0)
anova.gam(m3.DOC.T0)
gam.check(m3.DOC.T0, rep=1000)
draw(m3.DOC.T0)
concrvity(m3.DOC.T0)
par(mfrow = c(2, 2))
plot(m3.DOC.T0, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T0<-plot_difference(
  m1.DOC.T0,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T0.mod.plot<-
  plot_smooths(
  model = m3.DOC.T0,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-0") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting


######## T1 model
m1.DOC.T1<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.DOC.T1<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.DOC.T1<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.DOC.AIC<-AIC(m1.DOC.T1, m2.DOC.T1, m3.DOC.T1)
# best is smooth by factor

summary(m1.DOC.T1)
anova.gam(m1.DOC.T1)
gam.check(m1.DOC.T1, rep=1000)
draw(m1.DOC.T1)
concrvity(m1.DOC.T1)
par(mfrow = c(2, 2))
plot(m1.DOC.T1, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T1<-plot_difference(
  m1.DOC.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T1.mod.plot<-
  plot_smooths(
  model = m1.DOC.T1,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-10") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting


# effect of treatment, smoothing significant across both treatments
# DOC higher in  unburned, relative to burned


########## T2
m1.DOC.T2<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.DOC.T2<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.DOC.T2<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.DOC.AIC<-AIC(m1.DOC.T2, m2.DOC.T2, m3.DOC.T2)

# best is smooth by factor

summary(m1.DOC.T2)
anova.gam(m1.DOC.T2)
gam.check(m1.DOC.T2, rep=1000)
draw(m1.DOC.T2)
concrvity(m1.DOC.T2)
par(mfrow = c(2, 2))
plot(m1.DOC.T2, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T2<-plot_difference(
  m1.DOC.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T2.mod.plot<-
  plot_smooths(
  model = m1.DOC.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-31") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# NO effect of treatment, smoothing significant across both treatments
# DOC equivalent in burned and unburned 
# DOC more variable/wonky across gradient in burned


########## T3
m1.DOC.T3<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.DOC.T3<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.DOC.T3<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.DOC.AIC<-AIC(m1.DOC.T3, m2.DOC.T3, m3.DOC.T3)
# best by factor smooth

summary(m1.DOC.T3)
anova.gam(m1.DOC.T3)
gam.check(m1.DOC.T3, rep=1000)
draw(m1.DOC.T3)
concrvity(m1.DOC.T3)
par(mfrow = c(2, 2))
plot(m1.DOC.T3, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T3<-plot_difference(
  m1.DOC.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T3.mod.plot<-
  plot_smooths(
  model = m1.DOC.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-59") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# effect of treatment, smoothing significant across both treatments
# DOC higher in  burned vs. unburned


########## T4
m1.DOC.T4<-gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.DOC.T4<-gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.DOC.T4<-gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.DOC.AIC<-AIC(m1.DOC.T4, m2.DOC.T4, m3.DOC.T4)

# best is global
summary(m3.DOC.T4)
anova.gam(m3.DOC.T4)
gam.check(m3.DOC.T4, rep=1000)
draw(m3.DOC.T4)
concrvity(m3.DOC.T4)
par(mfrow = c(2, 2))
plot(m3.DOC.T4, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T4<-plot_difference(
  m1.DOC.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
DOC.T4.mod.plot<-
  plot_smooths(
  model = m3.DOC.T4,
  series = plant.mass..g
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-89") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# no effect of treatment, smoothing significant across both treatments
# DOC equivalent in  burned and unburned


mod.rep<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=5)

mod.DOC.df<- data.frame(mod.rep)

AIC.DOC<-bind_rows(T0.DOC.AIC, T1.DOC.AIC, T2.DOC.AIC, T3.DOC.AIC, T4.DOC.AIC)
AIC.DOC.mod<-cbind(mod.DOC.df, AIC.DOC)

write.csv(AIC.DOC.mod, "output/AIC models/AIC.DOC.csv")

DOC tables

Table: Results for DOC Time-0 (Day-0)

anova.gam(m3.DOC.T0)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                   edf Ref.df     F p-value
## s(plant.mass..g) 1.04   1.08 1.341    0.24

Table: Results for DOC Time-1 (Day-10)

anova.gam(m1.DOC.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.035   0.853
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   6.482  7.532 34.39  <2e-16
## s(plant.mass..g):Treatmentunburned 1.568  1.929 59.34  <2e-16

Table: Results for DOC Time-2 (Day-31)

anova.gam(m1.DOC.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.035   0.853
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   6.482  7.532 34.39  <2e-16
## s(plant.mass..g):Treatmentunburned 1.568  1.929 59.34  <2e-16

Table: Results for DOC Time-3 (Day-59)

anova.gam(m1.DOC.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 12.32 0.00182
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.051  2.532 94.00  <2e-16
## s(plant.mass..g):Treatmentunburned 2.202  2.714 56.55  <2e-16

Table: Results for DOC Time-4 (Day-89)

anova.gam(m3.DOC.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DOC..mg.L ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 1.928  2.385 29.8  <2e-16

DOC figures

Compile raw plots and model-diff plots for final figures, then plot the model difference.

###### compile the plots with effect plots

extract.legend <- get_legend(
  # create some space to the left of the legend
  DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


DOC.alltimes<-plot_grid(
  DOC.T0.mod.plot+ theme(legend.position = "none"),
  DOC.T1.mod.plot+ theme(legend.position = "none"),
  DOC.T2.mod.plot+ theme(legend.position = "none"),
  DOC.T3.mod.plot+ theme(legend.position = "none"),
  DOC.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/DOC.alltimes.mod.pdf", height=4, width=15)

DOC.alltimes
**Figure**. Dissolved organic carbon (DOC) concentration across time in treatments receiving burned and unburned plant material at the start of the experiment and four sampling periods after plant material added. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

Figure. Dissolved organic carbon (DOC) concentration across time in treatments receiving burned and unburned plant material at the start of the experiment and four sampling periods after plant material added. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

DOC.mod.diffs<-plot_grid(
  DOC.diff.T0+ theme(legend.position = "none")+ ggtitle("Day-0"),
  DOC.diff.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
  DOC.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  DOC.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  DOC.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8,8), ncol=5)
ggsave("figures/DOC.mod.diffs.pdf", height=4, width=14)

DOC.mod.diffs
**Figure**. Model effects from GAMs with differences between smoothers for DOC concentration across time in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Figure. Model effects from GAMs with differences between smoothers for DOC concentration across time in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

TDN

Using the above dataframes for DOC and TDN, analyze total N (total dissolved nitrogen (TDN)) and make plots in each timepoint, running models and making model-difference plots.

TN.df<-DOC.df

######## T0 model
m1.TN.T0<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T0", data = TN.df, method = "REML")

m2.TN.T0<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T0", data = TN.df, method = "REML")

m3.TN.T0<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T0", data = TN.df, method = "REML")

T0.TN.AIC<-AIC(m1.TN.T0, m2.TN.T0, m3.TN.T0)
# best smooth by factor

summary(m1.TN.T0)
anova.gam(m1.TN.T0)
gam.check(m1.TN.T0, rep=1000)
draw(m1.TN.T0)
concrvity(m1.TN.T0)
par(mfrow = c(2, 2))
plot(m1.TN.T0, all.terms = TRUE, page=1)

# model predictions
TN.diff.T0<-plot_difference(
  m1.TN.T0,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T0.mod.plot<-
  plot_smooths(
  model = m1.TN.T0,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-0") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# no treatment effect, compare to simplified model (p=0.901)

######## T1 model
m1.TN.T1<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T1", data = TN.df, method = "REML")

m2.TN.T1<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T1", data = TN.df, method = "REML")

m3.TN.T1<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T1", data = TN.df, method = "REML")

T1.TN.AIC<-AIC(m1.TN.T1, m2.TN.T1, m3.TN.T1)
#best global only

summary(m3.TN.T1)
anova.gam(m3.TN.T1)
gam.check(m3.TN.T1, rep=1000)
draw(m3.TN.T1)
concrvity(m3.TN.T1)
par(mfrow = c(2, 2))
plot(m3.TN.T1, all.terms = TRUE, page=1)

# model predictions
TN.diff.T1<-plot_difference(
  m1.TN.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T1.mod.plot<-
  plot_smooths(
  model = m3.TN.T1,
  series = plant.mass..g
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-10") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# TN smoother significant for unburned but not burned (p=0.007)


######## T2 model
m1.TN.T2<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T2", data = TN.df, method = "REML")

m2.TN.T2<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T2", data = TN.df, method = "REML")

m3.TN.T2<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T2", data = TN.df, method = "REML")

T2.TN.AIC<-AIC(m1.TN.T2, m2.TN.T2, m3.TN.T2)
#best global with treatment term

summary(m2.TN.T2)
anova.gam(m2.TN.T2)
gam.check(m2.TN.T2, rep=1000)
draw(m2.TN.T2)
concrvity(m2.TN.T2)
par(mfrow = c(2, 2))
plot(m2.TN.T2, all.terms = TRUE, page=1)

# model predictions
TN.diff.T2<-plot_difference(
  m1.TN.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T2.mod.plot<-
  plot_smooths(
  model = m2.TN.T2,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-31") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# Near treatment effect p=0.053, higher TN in unburned
# smoother signif: 0.042


######## T3 model
m1.TN.T3<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T3", data = TN.df, method = "REML")

m2.TN.T3<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T3", data = TN.df, method = "REML")

m3.TN.T3<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T3", data = TN.df, method = "REML")

T3.TN.AIC<-AIC(m1.TN.T3, m2.TN.T3, m3.TN.T3)
#best with smooth by factor term

summary(m1.TN.T3)
anova.gam(m1.TN.T3)
gam.check(m1.TN.T3, rep=1000)
draw(m1.TN.T3)
concrvity(m1.TN.T3)
par(mfrow = c(2, 2))
plot(m1.TN.T3, all.terms = TRUE, page=1)

# model predictions
TN.diff.T3<-plot_difference(
  m1.TN.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T3.mod.plot<-
  plot_smooths(
  model = m1.TN.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-59") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting

# Near treatment effect p=0.075, trend for higher TN in burned
# smoother signif: at <0.001 for both


######## T4 model
m1.TN.T4<-gam(TN..mg.L ~ Treatment + s(plant.mass..g, by=Treatment), 
              subset = Time.point=="T4", data = TN.df, method = "REML")

m2.TN.T4<-gam(TN..mg.L ~ Treatment + s(plant.mass..g), 
              subset = Time.point=="T4", data = TN.df, method = "REML")

m3.TN.T4<-gam(TN..mg.L ~ s(plant.mass..g), 
              subset = Time.point=="T4", data = TN.df, method = "REML")

T4.TN.AIC<-AIC(m1.TN.T4, m2.TN.T4, m3.TN.T4)
#best  with smooth by factor term

summary(m1.TN.T4)
anova.gam(m1.TN.T4)
gam.check(m1.TN.T4, rep=1000)
draw(m1.TN.T4)
concrvity(m1.TN.T4)
par(mfrow = c(2, 2))
plot(m1.TN.T4, all.terms = TRUE, page=1)

# model predictions
TN.diff.T4<-plot_difference(
  m1.TN.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## plot for the model output on rawdata
TN.T4.mod.plot<-
  plot_smooths(
  model = m1.TN.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=TN.df[(TN.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=TN..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +  
  coord_cartesian(ylim=c(0, 2)) +
  ggtitle("Day-89") +
  ylab("TN (mg/L)")  +
  xlab("plant material (g)") +
  Fig.formatting


# effect of treatment (higher TN in the burned) (p=0.020)
# significant smoother effect for burned treatment only (p=0.032)

mod.rep<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=5)

mod.TN.df<- data.frame(mod.rep)

AIC.TN<-bind_rows(T0.TN.AIC, T1.TN.AIC, T2.TN.AIC, T3.TN.AIC, T4.TN.AIC)
AIC.TN.mod<-cbind(mod.TN.df, AIC.TN)
write.csv(AIC.TN.mod, "output/AIC models/AIC.TN.mod.csv")

TDN tables

Results for TDN Time-0 (Day-0)

anova.gam(m1.TN.T0)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.879   0.357
## 
## Approximate significance of smooth terms:
##                                    edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned     1      1 0.009  0.9266
## s(plant.mass..g):Treatmentunburned   1      1 6.303  0.0186

Table: Results for TDN Time-1 (Day-10)

anova.gam(m3.TN.T1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 2.848  3.492 5.72 0.00274

Table: Results for TDN Time-2 (Day-31)

anova.gam(m2.TN.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 4.122  0.0532
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 3.207  3.921 4.87 0.00425

Table: Results for TDN Time-3 (Day-59)

anova.gam(m1.TN.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df   F p-value
## Treatment  1 3.5  0.0752
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value
## s(plant.mass..g):Treatmentburned   4.359  5.269 23.03  < 2e-16
## s(plant.mass..g):Treatmentunburned 2.457  3.022 10.52 0.000194

Table: Results for TDN Time-4 (Day-89)

anova.gam(m1.TN.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TN..mg.L ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 6.231  0.0196
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.417  2.973 3.613  0.0318
## s(plant.mass..g):Treatmentunburned 1.000  1.000 0.531  0.4731

TDN figures

Compile raw plots and model-diff plots for final figures, then plot the model difference.

###### compile the plots with effect plots

TN.mod.alltimes<-plot_grid(
  TN.T0.mod.plot+ theme(legend.position = "none"), 
  TN.T1.mod.plot+ theme(legend.position = "none"),
  TN.T2.mod.plot+ theme(legend.position = "none"),
  TN.T3.mod.plot+ theme(legend.position = "none"),
  TN.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/TN.mods.plots.pdf", height=4, width=13)

TN.mod.alltimes
**Figure**. Total dissolved nitrogen (TDN) concentration across time in treatments receiving burned and unburned plant material.  Lines represent best-fit GAMs with 95% confidence intervals.  Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

Figure. Total dissolved nitrogen (TDN) concentration across time in treatments receiving burned and unburned plant material. Lines represent best-fit GAMs with 95% confidence intervals. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

TN.mod.diffs<-plot_grid(
  TN.diff.T0+ theme(legend.position = "none")+ ggtitle("Day-0"),
  TN.diff.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
  TN.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  TN.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  TN.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/TN.mod.diffs.pdf", height=4, width=13)

TN.mod.diffs
**Figure**. Model effects from GAMs with differences between smoothers for TN across time in treatments receiving burned and unburned plant material. Model differences between the two factor-smoothers. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Figure. Model effects from GAMs with differences between smoothers for TN across time in treatments receiving burned and unburned plant material. Model differences between the two factor-smoothers. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

H2O Phosphorus

The phosphorus concentration (total dissolved phosphorus (TDP)) in water was only run at Time-2 (Day-31). Water was collected from each tank, and filtered using a GF/F.

  • We will make a plot of the raw data and run the model fitting, with smoother contrasts.
phos<-read.csv("data/Pyro_water.phosph.csv")

fac2<-c("Time.Point", "Treatment", "Tank") # make factors
phos[fac2]<-lapply(phos[fac2],factor)

phos<-na.omit(phos)

############# phosphorus T2: significant smoothers, not treatment overall
m1.T2.phos<- gam(TP.umol..l ~ Treatment +
            s(plant.mass..g, by=Treatment), data=phos, method="REML", family="gaussian")

m2.T2.phos<- gam(TP.umol..l ~ Treatment +
            s(plant.mass..g), data=phos, method="REML", family="gaussian")

m3.T2.phos<- gam(TP.umol..l ~
            s(plant.mass..g), data=phos, method="REML", family="gaussian")

AIC.T2.phos<-AIC(m1.T2.phos,m2.T2.phos, m3.T2.phos)
# factor smooth best

summary(m1.T2.phos)
anova.gam(m1.T2.phos)
gam.check(m1.T2.phos, rep=1000)
draw(m1.T2.phos)
concrvity(m1.T2.phos)
par(mfrow = c(1, 2))
plot(m1.T2.phos, all.terms = TRUE, page=1)

# model predictions
phos.diff.T2<-plot_difference(
  m1.T2.phos,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

## AIC table
mod.phos<-rep(c("Treatment + s(plant.mass..g, by=Treatment)", 
              "Treatment + s(plant.mass..g)", 
              "s(plant.mass..g)"), times=1)

mod.phos.df<- data.frame(mod.phos)

AIC.phos.mod<-cbind(mod.phos.df, AIC.T2.phos)
write.csv(AIC.phos.mod, "output/AIC models/AIC.phos.mod.csv")

Table: Results for total phosphorous across burned/unburned treatments at Day-31.

anova.gam(m1.T2.phos)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## TP.umol..l ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 1.329   0.267
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.924  3.525 154.7  <2e-16
## s(plant.mass..g):Treatmentunburned 2.930  3.371 124.6  <2e-16


###########  
#plot for the model output on rawdata
phos.T2.mod.plot<-
  plot_smooths(
  model = m1.T2.phos,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=phos, 
             aes(x=plant.mass..g, y=TP.umol..l, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("TP", ~(mu*mol/L), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


phos.plots<-plot_grid(phos.T2.mod.plot, phos.diff.T2, ncol=2,  rel_widths = c(8, 5), 
                      labels=c('A', 'B'), label_size=8)

phos.plots
**Figure**. (**A**) Total phosphorus concentration in water from burned and unburned treatments at Day-31, and (**B**) the difference between burned and unburned treatment smoothers.

Figure. (A) Total phosphorus concentration in water from burned and unburned treatments at Day-31, and (B) the difference between burned and unburned treatment smoothers.

ggsave("figures/phos.plots.pdf", height=5, width=8)

YSI

Import YSI data and produce plots of changes in O2% and net primary productivity (NPP) and respiration (Resp). The YSI data includes Temp, pH, dissolved oxygen (percent and concentration), and conductivity.
Here, we will pull in the raw data and make the new metrics NPP and Resp, determined from differences in DO% from dawn-dusk (NPP) and dusk-dawn (Resp) over a 24h period in each time period.

#load YSI data
YSI<-read.csv("data/Pyro_YSI.csv")

# fix date
YSI$Date<-as.character(YSI$Date)
YSI$Date<-as.POSIXct(YSI$Date, format="%m/%d/%y")
YSI$Date<-as.Date(YSI$Date, format="%m/%d/%Y")


####### Time 1 change in O2 ################

#separate time points
YSI.T1<- YSI[(YSI$Time.point=="T1"),]

#calculate NPP for T1
T1.Prod<-YSI.T1[(YSI.T1$Date == "2021-11-15"),] # dawn and dusk for 12h period
T1.Dawn1<-T1.Prod[(T1.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T1.Dusk<-T1.Prod[(T1.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T1.Dawn2<-YSI.T1[(YSI.T1$Date == "2021-11-16"),] # dawn-2 measurements, following AM

# make new dataframe
T1.O2<-(T1.Dawn1[,c(2,4:6)]) 
T1.O2$dawn1<-T1.Dawn1$DO.percent
T1.O2$dusk1<-T1.Dusk$DO.percent
T1.O2$dawn2<-T1.Dawn2$DO.percent

# R = dusk - dawn (PM to AM, O2 change of day 1)
# NPP = dusk - dawn (PM to AM, O2 change of day 2)

T1.O2<- mutate(T1.O2, 
                NPP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T1.O2<-T1.O2 %>% 
  arrange(Treatment, plant.mass..g) 


################ ################ ################
####### Time 2 change in O2 ################

#separate time points
YSI.T2<- YSI[(YSI$Time.point=="T2"),]

#calculate NPP for T2
T2.Prod<-YSI.T2[(YSI.T2$Date == "2021-12-06"),] # dawn and dusk for 12h period
T2.Dawn1<-T2.Prod[(T2.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T2.Dusk<-T2.Prod[(T2.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T2.Dawn2<-YSI.T2[(YSI.T2$Date == "2021-12-07"),] # dawn-2 measurements, following AM

# make new dataframe
T2.O2<-(T2.Dawn1[,c(2,4:6)]) 
T2.O2$dawn1<-T2.Dawn1$DO.percent
T2.O2$dusk1<-T2.Dusk$DO.percent
T2.O2$dawn2<-T2.Dawn2$DO.percent

T2.O2<- mutate(T2.O2, 
                NPP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T2.O2<-T2.O2 %>% 
  arrange(Treatment, plant.mass..g) 


################ ################ ################
####### Time 3 change in O2 ################

#separate time points
YSI.T3<- YSI[(YSI$Time.point=="T3"),]

#calculate NPP for T3
T3.Prod<-YSI.T3[(YSI.T3$Date == "2022-01-03"),] # dawn and dusk for 12h period
T3.Dawn1<-T3.Prod[(T3.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T3.Dusk<-T3.Prod[(T3.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T3.Dawn2<-YSI.T3[(YSI.T3$Date == "2022-01-04"),] # dawn-2 measurements, following AM

# make new dataframe
T3.O2<-(T3.Dawn1[,c(2,4:6)]) 
T3.O2$dawn1<-T3.Dawn1$DO.percent
T3.O2$dusk1<-T3.Dusk$DO.percent
T3.O2$dawn2<-T3.Dawn2$DO.percent

T3.O2<- mutate(T3.O2, 
                NPP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T3.O2<-T3.O2 %>% 
  arrange(Treatment, plant.mass..g) 

################ ################ ################
####### Time 3 change in O2 ################

#separate time points
YSI.T4<- YSI[(YSI$Time.point=="T4"),]

#calculate NPP for T4
T4.Prod<-YSI.T4[(YSI.T4$Date == "2022-02-02"),] # dawn and dusk for 12h period
T4.Dawn1<-T4.Prod[(T4.Prod$Dawn..Dusk == "dawn"),] # dawn-1 measurements
T4.Dusk<-T4.Prod[(T4.Prod$Dawn..Dusk == "dusk"),] # dusk measurements

T4.Dawn2<-YSI.T4[(YSI.T4$Date == "2022-02-03"),] # dawn-2 measurements, following AM

# make new dataframe
T4.O2<-(T4.Dawn1[,c(2,4:6)]) 
T4.O2$dawn1<-T4.Dawn1$DO.percent
T4.O2$dusk1<-T4.Dusk$DO.percent
T4.O2$dawn2<-T4.Dawn2$DO.percent

T4.O2<- mutate(T4.O2, 
                NPP=dusk1 - dawn1,
                R=dawn2 - dusk1) 

#sort
T4.O2<-T4.O2 %>% 
  arrange(Treatment, plant.mass..g) 

################ ################ ################
# combine T1  T2  T3 T4 timepoints
################ ################ ################
O2.tanks<-rbind(T1.O2,T2.O2, T3.O2, T4.O2)

cols<-c("Time.point", "Treatment", "Tank") # columns to make factors
O2.tanks[cols] <- lapply(O2.tanks[cols], factor) # make all these factors
O2.tanks$plant.mass..g<-as.numeric(O2.tanks$plant.mass..g)

DO and O2%

First, we will also run model fitting on the raw DO% data to apply the same approach for visualizing changes in oxygen across dawn-dusk-dawn measurements. We will then combine all these plots into multi-panel figures for NPP and Resp, and DO% for each point of measurement.

#TIME POINT 1**: Change in O~2~% from dissolved oxygen
#########################################################
##################################################################
# total oxygen % plot for the 3 time points (dawn-dusk-dawn)

m1.dawn1.T1<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.dawn1.T1<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.dawn1.T1<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.dawn1.AIC<-AIC(m1.dawn1.T1, m2.dawn1.T1, m3.dawn1.T1)
# global with treatment best


summary(m2.dawn1.T1)
anova.gam(m2.dawn1.T1)
gam.check(m2.dawn1.T1, rep=1000)
draw(m2.dawn1.T1)
concrvity(m2.dawn1.T1)
par(mfrow = c(2, 2))
plot(m2.dawn1.T1, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T1<-plot_difference(
  m1.dawn1.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T1.mod.plot<-
  plot_smooths(
  model = m2.dawn1.T1,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T1.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# treatment (p=0.0279) and smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T1<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.dusk1.T1<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.dusk1.T1<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.dusk1.AIC<-AIC(m1.dusk1.T1, m2.dusk1.T1, m3.dusk1.T1)
# model with treatment and global smooth best

summary(m2.dusk1.T1)
anova.gam(m2.dusk1.T1)
gam.check(m2.dusk1.T1, rep=1000)
draw(m2.dusk1.T1)
concrvity(m2.dusk1.T1)
par(mfrow = c(2, 2))
plot(m2.dusk1.T1, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T1<-plot_difference(
  m1.dusk1.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T1.mod.plot<-
  plot_smooths(
  model = m2.dusk1.T1,
  series = plant.mass..g,
  comparison= Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T1.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both treatments 


####### #### Dawn 2
m1.dawn2.T1<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.dawn2.T1<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.dawn2.T1<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.dawn2.AIC<-AIC(m1.dawn2.T1, m2.dawn2.T1, m3.dawn2.T1)
# treatment and global smooth best

summary(m2.dawn2.T1)
anova.gam(m2.dawn2.T1)
gam.check(m2.dawn2.T1, rep=1000)
draw(m2.dawn2.T1)
concrvity(m2.dawn2.T1)
par(mfrow = c(2, 2))
plot(m2.dawn2.T1, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T1<-plot_difference(
  m1.dawn2.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T1.mod.plot<-
  plot_smooths(
  model = m2.dawn2.T1,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  
  ggtitle("T1.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both (p<0.001)

#### group plots
O2.T1<-plot_grid(
  dawn1.T1.mod.plot+ theme(legend.position = "none"), 
  dusk1.T1.mod.plot+ theme(legend.position = "none"),
  dawn2.T1.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)
#TIME POINT 2**: Change in O~2~% from dissolved oxygen
#*############################################################
##############################################################################
# total oxygen % plot for the 3 time points (dawn-dusk-dawn)

m1.dawn1.T2<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.dawn1.T2<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.dawn1.T2<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.dawn1.AIC<-AIC(m1.dawn1.T2, m2.dawn1.T2, m3.dawn1.T2)
# factor by smooth best


summary(m1.dawn1.T2)
anova.gam(m1.dawn1.T2)
gam.check(m1.dawn1.T2, rep=1000)
draw(m1.dawn1.T2)
concrvity(m1.dawn1.T2)
par(mfrow = c(2, 2))
plot(m1.dawn1.T2, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T2<-plot_difference(
  m1.dawn1.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T2.mod.plot<-
  plot_smooths(
  model = m1.dawn1.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T2.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T2<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.dusk1.T2<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.dusk1.T2<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.dusk1.AIC<-AIC(m1.dusk1.T2, m2.dusk1.T2, m3.dusk1.T2)
# model with smooth by factor best

summary(m1.dusk1.T2)
anova.gam(m1.dusk1.T2)
gam.check(m1.dusk1.T2, rep=1000)
draw(m1.dusk1.T2)
concrvity(m1.dusk1.T2)
par(mfrow = c(2, 2))
plot(m1.dusk1.T2, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T2<-plot_difference(
  m1.dusk1.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T2.mod.plot<-
  plot_smooths(
  model = m1.dusk1.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T2.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both treatments 


####### #### Dawn 2
m1.dawn2.T2<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.dawn2.T2<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.dawn2.T2<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.dawn2.AIC<-AIC(m1.dawn2.T2, m2.dawn2.T2, m3.dawn2.T2)
# smooth by factor best

summary(m1.dawn2.T2)
anova.gam(m1.dawn2.T2)
gam.check(m1.dawn2.T2, rep=1000)
draw(m1.dawn2.T2)
concrvity(m1.dawn2.T2)
par(mfrow = c(2, 2))
plot(m1.dawn2.T2, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T2<-plot_difference(
  m1.dawn2.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T2.mod.plot<-
  plot_smooths(
  model = m1.dawn2.T2,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T2.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# smoother significant for both (p<0.001)

#### group plots
O2.T2<-plot_grid(
  dawn1.T2.mod.plot+ theme(legend.position = "none"), 
  dusk1.T2.mod.plot+ theme(legend.position = "none"),
  dawn2.T2.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)
#TIME POINT 3**: Change in O~2~% from dissolved oxygen
#############################################################
##############################################################################
#### Dawn1
m1.dawn1.T3<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.dawn1.T3<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.dawn1.T3<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.dawn1.AIC<-AIC(m1.dawn1.T3, m2.dawn1.T3, m3.dawn1.T3)
# factor by smooth best


summary(m1.dawn1.T3)
anova.gam(m1.dawn1.T3)
gam.check(m1.dawn1.T3, rep=1000)
draw(m1.dawn1.T3)
concrvity(m1.dawn1.T3)
par(mfrow = c(2, 2))
plot(m1.dawn1.T3, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T3<-plot_difference(
  m1.dawn1.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T3.mod.plot<-
  plot_smooths(
  model = m1.dawn1.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T3.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T3<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.dusk1.T3<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.dusk1.T3<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.dusk1.AIC<-AIC(m1.dusk1.T3, m2.dusk1.T3, m3.dusk1.T3)
# model with smooth by factor best

summary(m1.dusk1.T3)
anova.gam(m1.dusk1.T3)
gam.check(m1.dusk1.T3, rep=1000)
draw(m1.dusk1.T3)
concrvity(m1.dusk1.T3)
par(mfrow = c(2, 2))
plot(m1.dusk1.T3, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T3<-plot_difference(
  m1.dusk1.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T3.mod.plot<-
  plot_smooths(
  model = m1.dusk1.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T3.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect
# smoother significant for burned 


####### #### Dawn 2
m1.dawn2.T3<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.dawn2.T3<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.dawn2.T3<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.dawn2.AIC<-AIC(m1.dawn2.T3, m2.dawn2.T3, m3.dawn2.T3)
# smooth by factor best

summary(m1.dawn2.T3)
anova.gam(m1.dawn2.T3)
gam.check(m1.dawn2.T3, rep=1000)
draw(m1.dawn2.T3)
concrvity(m1.dawn2.T3)
par(mfrow = c(2, 2))
plot(m1.dawn2.T3, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T3<-plot_difference(
  m1.dawn2.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T3.mod.plot<-
  plot_smooths(
  model = m1.dawn2.T3,
  series = plant.mass..g,
  comparison=Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T3.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# global smoother significant(p<0.001)

#### group plots
O2.T3<-plot_grid(
  dawn1.T3.mod.plot+ theme(legend.position = "none"), 
  dusk1.T3.mod.plot+ theme(legend.position = "none"),
  dawn2.T3.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)
#TIME POINT 4: Change in O~2~% from dissolved oxygen
#
############################################################
##############################################################################
# total oxygen % plot for the 3 time points (dawn-dusk-dawn)

#### Dawn1
m1.dawn1.T4<-gam(dawn1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.dawn1.T4<-gam(dawn1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.dawn1.T4<-gam(dawn1 ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.dawn1.AIC<-AIC(m1.dawn1.T4, m2.dawn1.T4, m3.dawn1.T4)
# model with global best


summary(m3.dawn1.T4)
anova.gam(m3.dawn1.T4)
gam.check(m3.dawn1.T4, rep=1000)
draw(m3.dawn1.T4)
concrvity(m3.dawn1.T4)
par(mfrow = c(2, 2))
plot(m3.dawn1.T4, all.terms = TRUE, page=1)

# model predictions
dawn1.diff.T4<-plot_difference(
  m1.dawn1.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dawn1.T4.mod.plot<-
  plot_smooths(
  model = m3.dawn1.T4,
  series = plant.mass..g,
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=dawn1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T4.dawn1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting
  
# smoothers significant (p<0.001)

  
####### #### Dusk 1 
m1.dusk1.T4<-gam(dusk1 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.dusk1.T4<-gam(dusk1 ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.dusk1.T4<-gam(dusk1 ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.dusk1.AIC<-AIC(m1.dusk1.T4, m2.dusk1.T4, m3.dusk1.T4)
# model with smooth by factor best

summary(m1.dusk1.T4)
anova.gam(m1.dusk1.T4)
gam.check(m1.dusk1.T4, rep=1000)
draw(m1.dusk1.T4)
concrvity(m1.dusk1.T4)
par(mfrow = c(2, 2))
plot(m1.dusk1.T4, all.terms = TRUE, page=1)

# model predictions
dusk1.diff.T4<-plot_difference(
  m1.dusk1.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
dusk1.T4.mod.plot<-
  plot_smooths(
  model = m1.dusk1.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=dusk1, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T4.dusk1")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect
# smoother significant for burned 


####### #### Dawn 2
m1.dawn2.T4<-gam(dawn2 ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.dawn2.T4<-gam(dawn2 ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.dawn2.T4<-gam(dawn2 ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.dawn2.AIC<-AIC(m1.dawn2.T4, m2.dawn2.T4, m3.dawn2.T4)
# global smooth best

summary(m3.dawn2.T4)
anova.gam(m3.dawn2.T4)
gam.check(m3.dawn2.T4, rep=1000)
draw(m3.dawn2.T4)
concrvity(m3.dawn2.T4)
par(mfrow = c(2, 2))
plot(m3.dawn2.T4, all.terms = TRUE, page=1)

# model predictions
dawn2.diff.T4<-plot_difference(
  m1.dawn2.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned")),
)

###########  
#plot for the model output on rawdata
dawn2.T4.mod.plot<-
  plot_smooths(
  model = m3.dawn2.T4,
  series = plant.mass..g,
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=dawn2, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ggtitle("T4.dawn2")+
  coord_cartesian(ylim=c(-30, 150)) +
  ylab(expression(paste("O"[2],"%"))) +
  xlab("plant material (g)") +
  theme(legend.position = "right") +
  Fig.formatting

# global smoother significant(p<0.001)

#### group plots
O2.T4<-plot_grid(
  dawn1.T4.mod.plot+ theme(legend.position = "none"), 
  dusk1.T4.mod.plot+ theme(legend.position = "none"),
  dawn2.T4.mod.plot+ theme(legend.position = "none"),
  extract.legend, 
  rel_widths = c(8,8,8,3), ncol=4)

Combine and export all the O2 data with plot-difference and model AIC tables


#### model differences
O2.mod.diffs<-plot_grid(
  dawn1.diff.T1+ theme(legend.position = "none")+ ggtitle("T1-Dawn1"),
  dusk1.diff.T1+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T1+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  
  dawn1.diff.T2+ theme(legend.position = "none")+ ggtitle("T2-Dawn1"),
  dusk1.diff.T2+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T2+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  
  dawn1.diff.T3+ theme(legend.position = "none")+ ggtitle("T3-Dawn1"),
  dusk1.diff.T3+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T3+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  
  dawn1.diff.T4+ theme(legend.position = "none")+ ggtitle("T4-Dawn1"),
  dusk1.diff.T4+ theme(legend.position = "none")+ ggtitle("Dusk1"),
  dawn2.diff.T4+ theme(legend.position = "none")+ ggtitle("Dawn2"),
  rel_widths = c(8,8,8), ncol=3, nrow=4)

ggsave("figures/O2.mod.diffs.pdf", height=10, width=7)

#### model and raw data
O2.mods<-plot_grid(
  O2.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
  O2.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  O2.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  O2.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=1, nrow=4)

ggsave("figures/O2.mod.pdf", height=11, width=9)
O2.mods

#bind the AIC tables
AIC.O2<-bind_rows(T1.dawn1.AIC, T1.dusk1.AIC, T1.dawn2.AIC,
                  T2.dawn1.AIC, T2.dusk1.AIC, T2.dawn2.AIC,
                  T3.dawn1.AIC, T3.dusk1.AIC, T3.dawn2.AIC,
                  T4.dawn1.AIC, T4.dusk1.AIC, T4.dawn2.AIC)

# make a model column
mod.rep12<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=12)

mod.O2.df<- data.frame(mod.rep12)
#bind table

AIC.O2.mod<-cbind(mod.O2.df, AIC.O2)

write.csv(AIC.O2.mod, "output/AIC models/AIC.O2.mod.csv")
**Figure.** Changes in dissolved oxygen concentration (%) at dawn and dusk across the four experimental periods. Data here were used to calculate net ecosystem production and respiration. Lines represent best-fit GAMs with treatment-level 95% confidence intervals.  Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

Figure. Changes in dissolved oxygen concentration (%) at dawn and dusk across the four experimental periods. Data here were used to calculate net ecosystem production and respiration. Lines represent best-fit GAMs with treatment-level 95% confidence intervals. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

NPP and Resp

Generate dataframes for NPP and R change in O2. We will use NPP and Resp, run gam model fits, and produce individual figures for each time point.

  • First, we will use NPP models for productivity measurements.
####### Time 1

m1.NPP.T1<-gam(NPP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.NPP.T1<-gam(NPP ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.NPP.T1<-gam(NPP ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.NPP.AIC<-AIC(m1.NPP.T1, m2.NPP.T1, m3.NPP.T1)
# model with plot smooth by factor not different from reduced model, go with smooth by factor

summary(m2.NPP.T1)
anova.gam(m2.NPP.T1)
gam.check(m2.NPP.T1, rep=1000)
draw(m2.NPP.T1)
concrvity(m2.NPP.T1)
par(mfrow = c(2, 2))
plot(m2.NPP.T1, all.terms = TRUE, page=1)

mean.NPP.T1<-aggregate(NPP~Treatment, subset = Time.point=="T1", data = O2.tanks, FUN=mean)

#### see this https://cran.r-project.org/web/packages/tidymv/vignettes/plot-smooths.html
# The difference smooth is difference between the smooths of two conditions (two levels in a factor). 
# Portions of the difference smooth confidence interval that do not include 0 are shaded in red.

# model predictions
NPP.diff.T1<-plot_difference(
  m2.NPP.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
NPP.T1.mod.plot<-
  plot_smooths(
  model = m2.NPP.T1,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=NPP, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.110), smoothers significant (p<0.006)

####### Time 2

m1.NPP.T2<-gam(NPP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.NPP.T2<-gam(NPP ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.NPP.T2<-gam(NPP ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.NPP.AIC<-AIC(m1.NPP.T2, m2.NPP.T2, m3.NPP.T2)
# model with plot smooth by factor not different from reduced model, go with smooth by factor


summary(m2.NPP.T2)
anova.gam(m2.NPP.T2)
gam.check(m2.NPP.T2, rep=1000)
draw(m2.NPP.T2)
concrvity(m2.NPP.T2)
par(mfrow = c(2, 2))
plot(m2.NPP.T2, all.terms = TRUE, page=1)

mean.NPP.T2<-aggregate(NPP~Treatment, subset = Time.point=="T2", data = O2.tanks, FUN=mean)

## plot for the model output on rawdata
# model predictions
NPP.diff.T2<-plot_difference(
  m2.NPP.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

NPP.T2.mod.plot<-
  plot_smooths(
  model = m2.NPP.T2,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=NPP, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# treatment effect (p=0.007)
# smoother significant for burned (p=0.002) but not unburned (p=0.326)


####### Time 3
m1.NPP.T3<-gam(NPP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.NPP.T3<-gam(NPP ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.NPP.T3<-gam(NPP ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.NPP.AIC<-AIC(m1.NPP.T3, m2.NPP.T3, m3.NPP.T3)
# model with plot smooth by factor not different from reduced model, go with smooth by factor


summary(m2.NPP.T3)
anova.gam(m2.NPP.T3)
gam.check(m2.NPP.T3, rep=1000)
draw(m2.NPP.T3)
concrvity(m2.NPP.T3)
par(mfrow = c(2, 2))
plot(m2.NPP.T3, all.terms = TRUE, page=1)

mean.NPP.T3<-aggregate(NPP~Treatment, subset = Time.point=="T3", data = O2.tanks, FUN=mean)

# model predictions
NPP.diff.T3<-plot_difference(
  m2.NPP.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

#plot for the model output on rawdata
NPP.T3.mod.plot<-
  plot_smooths(
  model = m2.NPP.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=NPP, color=Treatment)) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# treatment effect (p=0.009)
# smoother significant for burned (p<0.001) but not unburned (p=0.053)


####### Time 4
m1.NPP.T4<-gam(NPP ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.NPP.T4<-gam(NPP ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.NPP.T4<-gam(NPP ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.NPP.AIC<-AIC(m1.NPP.T4, m2.NPP.T4, m3.NPP.T4)
# model with plot smooth by factor not different from reduced model, go with smooth by factor


summary(m1.NPP.T4)
anova.gam(m1.NPP.T4)
gam.check(m1.NPP.T4, rep=1000)
draw(m1.NPP.T4)
concrvity(m1.NPP.T4)
par(mfrow = c(2, 2))
plot(m1.NPP.T4, all.terms = TRUE, page=1)


# model predictions
NPP.diff.T4<-plot_difference(
  m1.NPP.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

#plot for the model output on rawdata
NPP.T4.mod.plot<-
  plot_smooths(
  model = m1.NPP.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=NPP, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-20, 50)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Net Production (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.118)
# smoother significant for burned (p=0.020) but not unburned (p=0.327)

mod.RNPP<-rep(c("~Treatment + s(plant.mass..g, by= Treatment)", 
              "~Treatment + s(plant.mass..g)", 
              "~s(plant.mass..g)"), times=4)

mod.RNPP.df<- data.frame(mod.RNPP)

AIC.NPP<-bind_rows(T1.NPP.AIC, T2.NPP.AIC, T3.NPP.AIC, T4.NPP.AIC)
AIC.NPP.mod<-cbind(mod.RNPP.df, AIC.NPP)
write.csv(AIC.NPP.mod, "output/AIC models/AIC.NPP.mod.csv")

Table: Results for NPP Time-1 (Day-10)

anova.gam(m1.NPP.T1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NPP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 2.756    0.11
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.478  3.046 13.44 2.3e-05
## s(plant.mass..g):Treatmentunburned 1.690  2.090  6.36 0.00585

Table: Results for NPP Time-2 (Day-31)

anova.gam(m1.NPP.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NPP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 8.479 0.00745
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.967  2.431 7.392 0.00227
## s(plant.mass..g):Treatmentunburned 1.000  1.000 1.002 0.32638

Table: Results for NPP Time-3 (Day-59)

anova.gam(m1.NPP.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NPP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 8.118 0.00948
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value
## s(plant.mass..g):Treatmentburned   3.414  4.166 8.236 0.000334
## s(plant.mass..g):Treatmentunburned 3.124  3.821 2.940 0.053008

Table: Results for NPP Time-4 (Day-89)

anova.gam(m1.NPP.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## NPP ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df    F p-value
## Treatment  1 2.62   0.118
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.757  3.382 3.717  0.0203
## s(plant.mass..g):Treatmentunburned 1.000  1.000 1.002  0.3268

Now, we will go through Respiration models and individual plots.

####### Time 1
m1.R.T1<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m2.R.T1<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

m3.R.T1<-gam(R ~ s(plant.mass..g), subset = Time.point=="T1", data = O2.tanks, method = "REML")

T1.R.AIC<-AIC(m1.R.T1, m2.R.T1, m3.R.T1)
# model with global best


summary(m3.R.T1)
anova.gam(m3.R.T1)
gam.check(m3.R.T1, rep=1000)
draw(m3.R.T1)
concrvity(m3.R.T1)
par(mfrow = c(2, 2))
plot(m3.R.T1, all.terms = TRUE, page=1)

# model predictions
R.diff.T1<-plot_difference(
  m1.R.T1,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T1.mod.plot<-
  plot_smooths(
  model = m3.R.T1,
  series = plant.mass..g,
) + 
  geom_point(data=T1.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.229), smoothers significant (p<0.001)

####### Time 2

m1.R.T2<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m2.R.T2<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

m3.R.T2<-gam(R ~ s(plant.mass..g), subset = Time.point=="T2", data = O2.tanks, method = "REML")

T2.R.AIC<-AIC(m1.R.T2, m2.R.T2, m3.R.T2)
# model with global + treatment best

summary(m2.R.T2)
anova.gam(m2.R.T2)
gam.check(m2.R.T2, rep=1000)
draw(m2.R.T2)
concrvity(m2.R.T2)
par(mfrow = c(2, 2))
plot(m2.R.T2, all.terms = TRUE, page=1)

mean.R.T2<-aggregate(R~Treatment, subset = Time.point=="T2", data = O2.tanks, FUN=mean)

# model predictions
R.diff.T2<-plot_difference(
  m2.R.T2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T2.mod.plot<-
  plot_smooths(
  model = m2.R.T2,
  series = plant.mass..g, 
  comparison= Treatment,
) + 
  geom_point(data=T2.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# slight treatment effect (p=0.085)
# smoother significant for both burned and unburned (p<0.008)


####### Time 3
m1.R.T3<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m2.R.T3<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

m3.R.T3<-gam(R ~ s(plant.mass..g), subset = Time.point=="T3", data = O2.tanks, method = "REML")

T3.R.AIC<-AIC(m1.R.T3, m2.R.T3, m3.R.T3)
# model smmoth by factor best


summary(m1.R.T3)
anova.gam(m1.R.T3)
gam.check(m1.R.T3, rep=1000)
draw(m1.R.T3)
concrvity(m1.R.T3)
par(mfrow = c(2, 2))
plot(m1.R.T3, all.terms = TRUE, page=1)

# model predictions
R.diff.T3<-plot_difference(
  m1.R.T3,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T3.mod.plot<-
  plot_smooths(
  model = m1.R.T3,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T3.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# treatment effect (p=0.027)
# smoother significant for burned and unnburned (p<0.001)


####### Time 4

m1.R.T4<-gam(R ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m2.R.T4<-gam(R ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

m3.R.T4<-gam(R ~ s(plant.mass..g), subset = Time.point=="T4", data = O2.tanks, method = "REML")

T4.R.AIC<-AIC(m1.R.T4, m2.R.T4, m3.R.T4)
# model with global + treatment best


summary(m1.R.T4)
anova.gam(m1.R.T4)
gam.check(m1.R.T4, rep=1000)
draw(m1.R.T4)
concrvity(m1.R.T4)
par(mfrow = c(2, 2))
plot(m1.R.T4, all.terms = TRUE, page=1)

# model predictions
R.diff.T4<-plot_difference(
  m1.R.T4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
R.T4.mod.plot<-
  plot_smooths(
  model = m1.R.T4,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=T4.O2, aes(x=plant.mass..g, y=R, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  coord_cartesian(ylim=c(-40, 10)) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab(expression(paste("Respiration (", Delta, "O"[2],"%)"))) +
  theme(legend.position = "right") +
  Fig.formatting

# no treatment effect (p=0.078)
# smoother significant for burned (p=0.004) but not unburned (p=0.965)


R.mod.plot<-plot_grid(
  R.T1.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-10"),
  R.T2.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-31"),
  R.T3.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-59"),
  R.T4.mod.plot+ theme(legend.position = "none")+ ggtitle("Day-89"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

#ggsave("figures/R.mod.plot.long.pdf", height=6, width=12)

#### model differences
R.mod.diffs<-plot_grid(
  R.diff.T1+ theme(legend.position = "none")+ ggtitle("R-Day-10"),
  R.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  R.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  R.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

#ggsave("figures/R.mod.diffs.pdf", height=3, width=10)

AIC.R<-bind_rows(T1.R.AIC, T2.R.AIC, T3.R.AIC, T4.R.AIC)
AIC.R.mod<-cbind(mod.RNPP.df, AIC.R)

write.csv(AIC.R.mod, "output/AIC models/AIC.R.mod.csv")

Table: Results for Resp Time-1 (Day-10)

anova.gam(m3.R.T1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F  p-value
## s(plant.mass..g) 2.520  3.097 12.8 2.32e-05

Table: Results for Resp Time-2 (Day-31)

anova.gam(m2.R.T2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 6.443  0.0186
## 
## Approximate significance of smooth terms:
##                    edf Ref.df  F  p-value
## s(plant.mass..g) 5.710  6.758 10 1.52e-05

Table: Results for Resp Time-3 (Day-39)

anova.gam(m1.R.T3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 5.669  0.0268
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F  p-value
## s(plant.mass..g):Treatmentburned   3.762  4.576 13.144 1.02e-05
## s(plant.mass..g):Treatmentunburned 3.274  4.000  7.775 0.000523

Table: Results for Resp Time-4 (Day-89)

anova.gam(m1.R.T4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## R ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df    F p-value
## Treatment  1 3.38  0.0784
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.927  3.587 5.293 0.00408
## s(plant.mass..g):Treatmentunburned 1.000  1.000 0.002 0.96531

Compile the NPP-Resp plots with model fits.

NPP.R.alltimes.long<-plot_grid(
  NPP.T1.mod.plot+ theme(legend.position = "none"),
  NPP.T2.mod.plot+ theme(legend.position = "none"),
  NPP.T3.mod.plot+ theme(legend.position = "none"),
  NPP.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  R.T1.mod.plot+ theme(legend.position = "none"),
  R.T2.mod.plot+ theme(legend.position = "none"),
  R.T3.mod.plot+ theme(legend.position = "none"),
  R.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3, 8,8,8,8,3), ncol=5, 
      labels=c('A','', '', '', '', 
               'B', '', '', '', ''), label_size=8)
ggsave("figures/NPP.R.alltimes.long.pdf", height=7, width=12)

NPP.R.alltimes.long
**Figure** (**A**) Net ecosystem productivity (NPP) and (**B**) respiration (Resp) in treatments receiving burned and unburned plant material across four sampling periods.  Lines represent best-fit generalized additive models (GAMs) with 95% confidence intervals.  Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

Figure (A) Net ecosystem productivity (NPP) and (B) respiration (Resp) in treatments receiving burned and unburned plant material across four sampling periods. Lines represent best-fit generalized additive models (GAMs) with 95% confidence intervals. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

NPP.R.mod.diffs<-plot_grid(
  NPP.diff.T1+ theme(legend.position = "none"),
  NPP.diff.T2+ theme(legend.position = "none"),
  NPP.diff.T3+ theme(legend.position = "none"),
  NPP.diff.T4+ theme(legend.position = "none"),
  R.diff.T1+ theme(legend.position = "none"),
  R.diff.T2+ theme(legend.position = "none"),
  R.diff.T3+ theme(legend.position = "none"),
  R.diff.T4+ theme(legend.position = "none"),
  rel_widths = c(8,8,8,8,8,8,8,8), ncol=4)
ggsave("figures/NPP.R.mod.diffs.alt.pdf", height=7, width=12)

NPP.R.mod.diffs
**Figure**. Model effects from GAMs with differences between smoothers for oxygen concentration (% O2) across four experimental time points in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Figure. Model effects from GAMs with differences between smoothers for oxygen concentration (% O2) across four experimental time points in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Isotopes

Isotopes and C:N (carbon to nitrogen molar concentrations) for starting materials used in the experiment and plankton fractions sampled at time 1 (Day-10) and time 2 (Day-31).


#########  Time 1 and Time 2
topes<-read.csv("data/Isotopes/Pyro_isotopes.csv")
topes$C.N <-(topes$Total.C..ug/12)/(topes$Total.N..ug/14) # C mol : N mol

cols<-c("Time.point", "Treatment", "Type", "Tank") # columns to make factors
topes[cols] <- lapply(topes[cols], factor) # make all these factors

##### make data frames
# treatment data df
topes.trt<-topes[(topes$Treatment=="burned" | topes$Treatment=="unburned"),]
topes.trt<-droplevels(topes.trt)
topes.trt$Type<-factor(topes.trt$Type, 
                         levels=c("plankton", "POM"))

# control and start plant materials df
topes.controls<-topes[!(topes$Treatment=="burned" | topes$Treatment=="unburned"),]

Starting materials

We will turn our attention to the starting materials, the burned and unburned 15N-labeled sage and non-labeled-willow (and non-labeled stock plankton) from the beginning of the experiment to see what their C:N values and isotope values are. We define these as “control” samples and include the tin blanks, plankton, and starting materials.

  • The code chunk below structures contrasts and runs non-parametric tests and linear models testing for differences between sample types.
  • Make a combined figure of nitrogen isotope values and C:N values.
########## ########## ########## 
# run some stats to see how the control material differs from each other

# make a burning treatment
topes.controls$Burn.Trt=as.factor(word(topes.controls$Type, -1, sep="[.]"))

# make a plant name
topes.controls$Plant=as.factor(word(topes.controls$Type, 1, sep="[.]"))

### test some models on controls
# remove blanks and plankton, keeping only plants
topes.controls.plants<-topes.controls[!(topes.controls$Plant=="blank" |
                                          topes.controls$Plant=="plankton" ),]
topes.controls.plants$Plant<-droplevels(topes.controls.plants$Plant)

# keep only non-enriched samples (remove sage)
topes.controls.non.enrich<-topes.controls[!(topes.controls$Plant=="blank" |
                                          topes.controls$Plant=="sage" ),]
topes.controls.non.enrich$Plant<-droplevels(topes.controls.non.enrich$Plant)


######## test plant species differences
mwu(topes.controls.plants, d15N, Plant)
# d15N sage and willow differ (p<0.001)
mwu(topes.controls.plants, C.N, Plant) 
# C.N sage and willow differ (p=0.013)

par(mfrow=c(1,2))
boxplot(d15N~Plant, data=topes.controls.plants)
boxplot(C.N~Plant, data=topes.controls.plants)


######## test difference between willow and plankton 
mwu(topes.controls.non.enrich, d15N, Plant)
 # d15N plankton and willow differ (p<0.001)
mwu(topes.controls.non.enrich, C.N, Plant) # C:N plankton and willow differ (p<0.001)
par(mfrow=c(1,2))
boxplot(d15N~Plant, data=topes.controls.non.enrich)
boxplot(C.N~Plant, data=topes.controls.non.enrich)


############# separate plant dfs
#### Sage d15N
topes.controls.sage<-topes.controls[(topes.controls$Plant=="sage"),]
topes.controls.sage$Plant<-droplevels(topes.controls.sage$Plant)
topes.controls.sage$Burn.Trt<-droplevels(topes.controls.sage$Burn.Trt)

# how do different types of sage compare across burn/unburn
# first, no difference between burned or very burned sage
anova(lm(d15N~Burn.Trt, data=topes.controls.sage[!(topes.controls.sage$Burn.Trt=="unburned"),]))

# convert to just 2 levels, no difference here either
topes.controls.sage$Burn.Unb<-ifelse(topes.controls.sage$Burn.Trt=="burned", "burned",
                                     ifelse(topes.controls.sage$Burn.Trt=="very burned", "burned",
                                     "unburned"))


mod.sage<-lm(d15N~Burn.Trt, data=topes.controls.sage) # keep at 3 levels
anova(mod.sage) 
# no difference in d15N for burned, unburned, very burned sage (p=0.423)

#### Sage C.N
mod.sage.CN<-lm(C.N~Burn.Trt, data=topes.controls.sage)
anova(mod.sage.CN) 
posthoc<-emmeans(mod.sage.CN, ~Burn.Trt)
multcomp::cld(posthoc, Letters=letters)
# Sage: difference in unburned, burned, very burned for C:N

par(mfrow=c(1,2))
boxplot(d15N~Burn.Trt, data=topes.controls.sage)
boxplot(C.N~Burn.Trt, data=topes.controls.sage)



###########################
#### Willow d15N
topes.controls.will<-topes.controls[(topes.controls$Plant=="willow"),]
topes.controls.will$Plant<-droplevels(topes.controls.will$Plant)
topes.controls.will$Burn.Trt<-droplevels(topes.controls.will$Burn.Trt)

mod<-lm(C.N~Burn.Trt, data=topes.controls.will)
anova(mod) 
# no difference in burned/unburned willow d15N

#### Willow C.N
mod<-lm(C.N~Burn.Trt, data=topes.controls.will)
anova(mod) 
# Willow: no difference in unburned, burned C:N (p=0.061)

par(mfrow=c(1,2))
boxplot(d15N~Burn.Trt, data=topes.controls.will)
boxplot(C.N~Burn.Trt, data=topes.controls.will)


######### make summary dfs
# summarize by plants
d15N.plant<-aggregate(d15N~Plant, topes.controls, FUN=mean)
d15N.plantSD<-aggregate(d15N~Plant, topes.controls, FUN=sd)
d15N.plant[3]<-d15N.plantSD[2]
colnames(d15N.plant)<-c("Plant", "d15N", "SD")

CN.plant<-aggregate(C.N~Plant, topes.controls, FUN=mean)
CN.plantSD<-aggregate(C.N~Plant, topes.controls, FUN=sd)
CN.plant[3]<-CN.plantSD[2]
colnames(CN.plant)<-c("Plant", "C.N", "SD")

# summary df d15N
d15N.cont<-aggregate(d15N~Type, topes.controls, FUN=mean)
d15N.cont.n<-aggregate(d15N~Type, topes.controls, FUN=length)
d15N.cont.SD<-aggregate(d15N~Type, topes.controls, FUN=sd)
d15N.cont[3]<- d15N.cont.SD[2]
d15N.cont[4]<- d15N.cont.n[2]
colnames(d15N.cont)<-c("Type", "d15N", "SD", "n")

# summary df control C:N
CN.cont<-aggregate(C.N~Type, topes.controls, FUN=mean)
CN.cont.n<-aggregate(C.N~Type, topes.controls, FUN=length)
CN.cont.SD<-aggregate(C.N~Type, topes.controls, FUN=sd)
CN.cont[3]<- CN.cont.SD[2]
CN.cont[4]<- CN.cont.n[2]
colnames(CN.cont)<-c("Type", "C.N", "SD", "n")
# make boxplots of control sample d15N and C:N
########## control plots
# set levels
topes.controls$Type<-factor(topes.controls$Type, 
                         levels=c("blank", "plankton.stock", 
                                  "willow.unburned", "willow.burned",
                                  "sage.unburned", "sage.burned", 
                                  "sage.veryburned", "sage.stem.burned"))

#### controls d15N boxplot
iso.plot.control.d15N<-ggplot(data=topes.controls,  aes(x=Type, y=d15N, fill=Type)) +
  geom_boxplot() +
  geom_point(pch = 21, position = position_jitterdodge(), alpha=0.6)+
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) +
  scale_fill_manual(values = c("azure2", "cornflowerblue",
                                "darkgoldenrod1", "indianred2",
                                "aquamarine3", "antiquewhite3",
                                "darkolivegreen4", "lightsalmon")) +
  xlab("control types") +  Fig.formatting +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 


###### control C:N
topes.controls.CN<-topes.controls %>% drop_na(C.N) # drop the NAs for C.N, makes plotting problematic

iso.plot.control.CN<-ggplot(data=topes.controls.CN,  aes(x=Type, y=C.N, fill=Type)) +
  geom_boxplot() +
  geom_point(pch = 21, position = position_jitterdodge(), alpha=0.6)+
  ylab("C:N") +
  scale_fill_manual(values = c("cornflowerblue",
                                "darkgoldenrod1", "indianred2",
                                "aquamarine3", "antiquewhite3",
                                "darkolivegreen4", "lightsalmon")) +
  xlab("control types") +  Fig.formatting +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


####### combine plots
# legend
extract.legend.cont <- get_legend(
  # create some space to the left of the legend
  iso.plot.control.d15N + theme(legend.box.margin = margin(0, 0, 0, 10)))


## combine 
control.iso.alltime<- 
  plot_grid(iso.plot.control.d15N + theme(legend.position = "none"),
            iso.plot.control.CN + theme(legend.position = "none"),
            extract.legend.cont, rel_widths = c(8,8,3), ncol=3, labels=c('A', 'B'), label_size=8)

control.iso.alltime
**Figure**. (**A**) Nitrogen isotope values and (**B**) C:N ratio for experimental controls (tin blanks), stock plankton, and burned or unburned plant material (willow, sage).

Figure. (A) Nitrogen isotope values and (B) C:N ratio for experimental controls (tin blanks), stock plankton, and burned or unburned plant material (willow, sage).

ggsave("figures/iso.controls.pdf", encod="MacRoman", height=5, width=10)

Plankton C:N

We will test C:N vs. plant material models and generate plot, fit with GAMs

  • Here we will test if the C:N was consistent across time, or what factors may influence C:N. This is important to verify assumptions.
  • We are using 15N transfer between trophic levels, assuming our measurements of efficiency in nitrogen transfer reflect carbon and energy transfer between trophic levels.

We will run through the model fitting and test for differences in the Treatments (Burned and Unburned plant material), Plankton Type (either > or < 63 μm), and Timepoint (Time 1 or Time 2).

############# all plankton T1 and T2
m1.CN <- gam(C.N ~ Treatment + Type + Time.point +
            s(plant.mass..g, by=Treatment),
            data=topes.trt, method="REML", family="gaussian")

m2.CN <- gam(C.N ~ Treatment + Type + Time.point +
            s(plant.mass..g), 
            data=topes.trt, method="REML", family="gaussian")

m3.CN <- gam(C.N ~
            s(plant.mass..g), 
            data=topes.trt, method="REML", family="gaussian")

m4.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g, by=Time.point),
            data=topes.trt, method="REML", family="gaussian")

m5.CN <- gam(C.N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment),
            data=topes.trt, method="REML", family="gaussian")

#best model here
m6.CN <- gam(C.N ~ Type +
            s(plant.mass..g, by=Treatment),
            data=topes.trt, method="REML", family="gaussian")

m7.CN <- gam(C.N ~ Treatment +
            s(plant.mass..g, by=Type),
            data=topes.trt, method="REML", family="gaussian")

m8.CN <- gam(C.N ~ Type +
            s(plant.mass..g, by=Type),
            data=topes.trt, method="REML", family="gaussian")


AIC.CN<-AIC(m1.CN, m2.CN, m3.CN, m4.CN, m5.CN, m6.CN, m7.CN, m8.CN)
## additive model best fit, but no treatment or type effect

summary(m6.CN)
anova.gam(m6.CN)
gam.check(m6.CN, rep=1000)
draw(m6.CN)
concrvity(m6.CN)
par(mfrow = c(1, 2))
plot(m6.CN, all.terms = TRUE, page=1)

# model for smoothing
msmooth.CN<-gam(C.N ~ Type +
            s(plant.mass..g, by=Type),
            data=topes.trt, method="REML", family="gaussian")

# model predictions
CN.diff<-plot_difference(
  m1.CN,
  series = plant.mass..g,
  difference = list(Type = c("plankton", "POM"))
)

#####

# linear model approach
CN.all.mod<-lm(C.N~ Treatment+Type+Time.point, data=topes.trt, na.action=na.exclude) 
print(Anova(CN.all.mod, type=2), digits=5)
posthoc<-emmeans(CN.all.mod, ~Type)
multcomp::cld(posthoc, Letters=letters)

##### plot for the model output on rawdata
CN.mod.plot.timepooled<-
  plot_smooths(
  model = msmooth.CN,
  series = plant.mass..g,
  comparison = Type)  + 
  theme(legend.position = "none") +
  geom_point(data=topes.trt, 
             aes(x=plant.mass..g, y=C.N, color=Type, fill=Type)) +
  scale_fill_manual(values = c("deepskyblue4", "darkseagreen"), guide='none') +
  scale_color_manual(name="Plankton", values = c("deepskyblue4", "darkseagreen"), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  theme(legend.position = "right") +
  ggtitle("Days 10 and 31") + 
  coord_cartesian(ylim=c(0, 20)) +
  ylab("C:N") +
  xlab("plant material (g)") +
  Fig.formatting 

ggsave("figures/CN.mod.plot.timepooled.pdf", height=4, width=5, encod="MacRoman")


############
### All time C.N boxplot
CNbox.all.time<-ggplot(topes.trt, aes(x=Treatment, y=C.N, fill=Type)) + 
  geom_boxplot() +
  geom_point(pch = 21, position = position_jitterdodge(), alpha=0.6) +
  scale_fill_manual(name="Plankton", values = c("deepskyblue4", "darkseagreen"), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  coord_cartesian(ylim=c(0, 20)) +
  ylab("C:N") +
  Fig.formatting

ggsave("figures/CNbox.all.time.pdf", height=4, width=5, encod="MacRoman")


C.Nboxplots<- plot_grid(CN.mod.plot.timepooled, 
                    CNbox.all.time,
                    rel_widths = c(8,8), ncol=2, labels=c('A', 'B'), label_size=8)
C.Nboxplots
**Figure**. (**A**) Plankton C:N along the plant material gradient pooled across days (10 and 31) and treatments (burned and unburned), and (**B**) plankton C:N in treatment tanks receiving burned and unburned plant material.

Figure. (A) Plankton C:N along the plant material gradient pooled across days (10 and 31) and treatments (burned and unburned), and (B) plankton C:N in treatment tanks receiving burned and unburned plant material.

Table: Results for C:N across size fractions pooled at Days-10 and 31.

anova.gam(m6.CN)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## C.N ~ Type + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##      df     F  p-value
## Type  1 22.58 6.04e-06
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value
## s(plant.mass..g):Treatmentburned   4.426  5.346 12.421 < 2e-16
## s(plant.mass..g):Treatmentunburned 1.650  2.039  6.618 0.00182

Table: Results for C:N across burned/unburned samples and size fractions pooled at Days-10 and 31.

print(Anova(CN.all.mod, type=2), digits=5)
## Anova Table (Type II tests)
## 
## Response: C.N
##            Sum Sq  Df F value    Pr(>F)    
## Treatment    0.08   1  0.0116 0.9145723    
## Type        98.50   1 13.7615 0.0003201 ***
## Time.point   0.55   1  0.0765 0.7825313    
## Residuals  830.26 116                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Trophic Transfer Efficiency

We will not use a two-member mixing model to calculate %sage-^15^N in the plankton. The %sage-15N is a proxy for trophic transfer efficiency, using the transfer of N from inorganic forms into phytoplankton and the incorporation of this N into higher trophic levels (zooplankton consumers). Here, the higher the enrichment of 15N in the plankton, the more the source of N reflects terrestrial material (i.e., 15N-labeled sage) and the more allochthonous nutrition is supplying N for plankton consumers.

  • Overall: Higher %sage-15N in plankton = more allochthonous/terrestrial subsidies N into plankton
  • Autochthonous end-member: we used the value of 11 per mil, which reflects the start values of plankton/POM/willow.
  • Terrestrial subsidy end-member: we used δ15N value of sage (no differed in burned vs. unburned) of 296 permil.
  • we recognize that the absolute contribution of terrestrial subsidies to plankton would be higher, since tanks received both sage and willow, but willow was not labeled.

To account for nitrogen addition, we will run the analysis with the mean %N for burned (1.296545 %N) and unburned (1.177917 %N) plant material x the mass of the detritus added to the tanks. This N-based x-axis represents the total N added to the tanks. To allow comparison to the raw data as δ15N and total plant material added, we also perform the same analyses and supply these figures in the supplement.

#### Format and run the mixing model for % sage.

# mixing model
head(topes.trt)
topes.trt<-droplevels(topes.trt)


### values for controls
d15N.cont # summary mean d15N by all controls
d15N.plant # summary by plants, # stock plankton 11, # sage 296, #  willow 13

# summary mean atom percent enrichment
F.cont<-aggregate(at.P..15N  ~ Type, topes.controls, mean)
F.plant<-aggregate(at.P..15N~Plant, topes.controls, FUN=mean)
# sage ~0.475
# willow 0.371

# 2 source mixing model (Post 2002), used d15N values here

# alpha = percent Sage from food web 1
# %Sage = (d15N sample - d15N base 2 [i.e., no-label food])/ (d15N sage food 1 - source 2)
#  d15N values of base 2 = 11 permil for algae/plankton stock
#  d15N value of base 1 = 296 permil for sage

# framed differently from Robinson 2001, TREE
# xtracer = frction of tracer
# Xtracer = (d15N-sample - d15N background) / (d15N-tracer - d15N-background)

topes.trt$percent.sage<-(topes.trt$d15N-12)/(296-11)*100

# unicode text for micrometer = \u03BC, use this in legend



########  15N-sage and N added
# summary from the elemental analysis
plant.nut<-read.csv("data/Pyro_plant material_elemental.csv")

cols<-c("type", "plant", "treatment") # columns to make factors
plant.nut[cols] <- lapply(plant.nut[cols], factor) # make all these factors

plant.sum.trt<-aggregate(N~treatment, plant.nut, FUN=mean)



#### looking at %N from elemental analysis
N.box.percent.plant.material<- ggplot(plant.nut, aes(x=treatment:plant, y=N, fill=treatment)) +
  geom_boxplot(alpha=0.7) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  geom_dotplot(binaxis='y', stackdir='center', alpha=0.5, dotsize=0.5,
               position=position_dodge(0.75))+
  xlab("Treatment")+
  ylab("%N")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

N.box.percent.plant.material
ggsave("figures/percent.N.box.trt.pdf", encod="MacRoman", height=4, width=4)
 

######### calculate %N from total plant biomass

# if all things equal, use the %N of sage and willow (stem and leaf) to determine the g of N added
topes.trt$plant.mass..g.N<- ifelse(topes.trt$Treatment =="burned",
                                   topes.trt$plant.mass..g*(1.296545/100),
                                   topes.trt$plant.mass..g*(1.177917/100))



######### plot relationship between N-grams and total plant biomass-grams
biomass.plant.N<-ggplot(topes.trt, aes(x=plant.mass..g, y=plant.mass..g.N, color=Treatment))+
  geom_point()+
  geom_line() + theme_classic()

biomass.plant.N
ggsave("figures/biomass.plant.N.pdf", encod="MacRoman", height=3, width=4)
############# all plankton T1

# test models T1
m1.T1.sage <- gam(percent.sage ~ Treatment + Type +
                    s(plant.mass..g.N, by=Treatment), 
                  subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m2.T1.sage <- gam(percent.sage ~ Treatment + Type +
                         s(plant.mass..g.N), 
                       subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m3.T1.sage <- gam(percent.sage ~
                    s(plant.mass..g.N), 
                  subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

AIC.sage.T1<-AIC(m1.T1.sage, m2.T1.sage, m3.T1.sage)

# anova for best model
anova.gam(m1.T1.sage)

# smooth fit for plot
msmooth.T1<- gam(percent.sage ~ Treatment +
                        s(plant.mass..g.N, by=Treatment), 
                      subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")
# model predictions
per.Sage.diff.T1<-plot_difference(
  msmooth.T1,
  series = plant.mass..g.N,
  difference = list(Treatment = c("burned", "unburned"))
)

#plot
per.Sage.T1.mod.plot<-
  plot_smooths(
    model = msmooth.T1,
    series = plant.mass..g.N,
    comparison = Treatment
  )  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T1"),], 
             aes(x=plant.mass..g.N, y=percent.sage, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("% Sage")+
  xlab("plant N added (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 100)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


### T2
# test models
m1.T2.sage <- gam(percent.sage ~ Treatment + Type +
                         s(plant.mass..g.N, by=Treatment), 
                       subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m2.T2.sage <- gam(percent.sage ~ Treatment + Type +
                         s(plant.mass..g.N), 
                       subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m3.T2.sage <- gam(percent.sage ~
                         s(plant.mass..g.N), 
                       subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

AIC.sage.T2<-AIC(m1.T2.sage, m2.T2.sage, m3.T2.sage)

# anova for best model
anova.gam(m1.T2.sage)


# smooth fit for plot
msmooth.T2<- gam(percent.sage ~ Treatment +
                   s(plant.mass..g.N, by=Treatment), 
                 subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

# model predictions
per.Sage.diff.T2<-plot_difference(
  msmooth.T2,
  series = plant.mass..g.N,
  difference = list(Treatment = c("burned", "unburned"))
)


# plot
per.Sage.T2.mod.plot<-
  plot_smooths(
    model = msmooth.T2,
    series = plant.mass..g.N,
    comparison = Treatment
  )  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T2"),], 
             aes(x=plant.mass..g.N, y=percent.sage, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab("% Sage")+
  xlab("plant N added (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 100)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


# effect of treatment (p<0.001) but not Type (p=0.321)
# smoother significant for both terms


### model plots for percent sage mixing model

## AIC table
mod.sag.top<-rep(c( "Treatment + Type + s(plant.mass..g.N, by=Treatment)", 
              "Treatment + Type + s(plant.mass..g.N)",
              "s(plant.mass..g.N)"), times=2)

mod.sag.df<- data.frame(mod.sag.top)
AIC.sag.topes<-bind_rows(AIC.sage.T1, AIC.sage.T2)

AIC.sag.mod<-cbind(mod.sag.df, AIC.sag.topes)
write.csv(AIC.sag.mod, "output/AIC models/AIC.sag.mod.csv")

Table: Results for %sage-15N across burned/unburned treatments at Day-10.

anova.gam(m1.T1.sage)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## percent.sage ~ Treatment + Type + s(plant.mass..g.N, by = Treatment)
## 
## Parametric Terms:
##           df      F  p-value
## Treatment  1 51.193 3.55e-09
## Type       1  1.704    0.198
## 
## Approximate significance of smooth terms:
##                                        edf Ref.df     F p-value
## s(plant.mass..g.N):Treatmentburned   3.537  4.301 136.7  <2e-16
## s(plant.mass..g.N):Treatmentunburned 3.791  4.584 178.3  <2e-16

Table: Results for %sage-15N across burned/unburned treatments at Day-31.

anova.gam(m1.T2.sage)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## percent.sage ~ Treatment + Type + s(plant.mass..g.N, by = Treatment)
## 
## Parametric Terms:
##           df      F  p-value
## Treatment  1 21.952 2.25e-05
## Type       1  1.003    0.322
## 
## Approximate significance of smooth terms:
##                                        edf Ref.df      F p-value
## s(plant.mass..g.N):Treatmentburned   4.631  5.542  80.65  <2e-16
## s(plant.mass..g.N):Treatmentunburned 3.301  4.021 127.03  <2e-16

Generate the combined model plots for %-sage-15N, and plot the model difference.

  • this plot is using the total plant N added (x-axis)
  • the total-plant N added to tanks = mean %N for burned or unburned detritus x total biomass added
# legend
extract.legend.mix <- get_legend(
  # create some space to the left of the legend
  per.Sage.T2.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


## combine 
sage.mix.model<- plot_grid(per.Sage.T1.mod.plot + theme(legend.position = "none"), 
                    per.Sage.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend.mix,
                    rel_widths = c(8,8,3), ncol=3, labels=c('A', 'B', ''), label_size=8)
sage.mix.model
**Figure**. Trophic transfer into particulate organic material (< 63 μm, circles) and zooplankton (> 63 μm, triangles) as the % sage-^15^N from a two-source mixing model as a metric for terrestrial subsidies at Days-10 and 31. Plant biomass nitrogen (x-axis) was calculated from the %N of and total biomass of plant burned and unburned plant detritus added to tanks receiving.  Lines represent best-fit generalized additive models (GAMs) with treatment-level 95% confidence intervals.

Figure. Trophic transfer into particulate organic material (< 63 μm, circles) and zooplankton (> 63 μm, triangles) as the % sage-15N from a two-source mixing model as a metric for terrestrial subsidies at Days-10 and 31. Plant biomass nitrogen (x-axis) was calculated from the %N of and total biomass of plant burned and unburned plant detritus added to tanks receiving. Lines represent best-fit generalized additive models (GAMs) with treatment-level 95% confidence intervals.


ggsave("figures/Isotope.mixmodel.pdf", encod="MacRoman", height=4, width=8)

# and plot difference
sage.mix.model.diff<- plot_grid(per.Sage.diff.T1 + theme(legend.position = "none"), 
                    per.Sage.diff.T2 + theme(legend.position = "none"),
                    rel_widths = c(8,8), ncol=2)
sage.mix.model.diff
**Figure** Model effects from GAMs with differences between smoothers for % sage-derived 15N at Day-10 and Day-31 in tanks receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Figure Model effects from GAMs with differences between smoothers for % sage-derived 15N at Day-10 and Day-31 in tanks receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

ggsave("figures/Isotope.mix.plotdiff.pdf", encod="MacRoman", height=4, width=8)

We also ran model fits and analysis on the raw isotope data for the plankton (δ15N) outside of the mixing model. See supplemental plots for output and code below.


## make d15N plots as well -- these follow the % sage, but are informative with the control plots to see the d15N of plankton and the 2 end members.

#### **** this plot is with total plant biomass and d15N **** ####
# other plot was %-sage-15N and plant-N-biomass

#### d15N isotope plots for treatments

############# all plankton T1
m1.T1.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m2.T1.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

m3.T1.d15N <- gam(d15N ~
            s(plant.mass..g), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

AIC.d15N.T1<-AIC(m1.T1.d15N, m2.T1.d15N, m3.T1.d15N)
## additive model best fit

summary(m1.T1.d15N)
anova.gam(m1.T1.d15N)
gam.check(m1.T1.d15N, rep=1000)
draw(m1.T1.d15N)
concrvity(m1.T1.d15N)
par(mfrow = c(1, 2))
plot(m1.T1.d15N, all.terms = TRUE, page=1)

# model predictions
d15N.diff.T1<-plot_difference(
  m1.T1.d15N,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
# model for smoothing
msmooth.T1.d15N<- gam(d15N ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T1", data=topes.trt, method="REML", family="gaussian")

#plot for the model output on rawdata
d15N.T1.mod.plot<-
  plot_smooths(
  model = msmooth.T1.d15N,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=d15N, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) +
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 250)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

# overall an effect of burning with both smoothers being significant by treatment
# no effect of type, POM and plankton with similar d15N values


############# all plankton T2
m1.T2.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m2.T2.d15N <- gam(d15N ~ Treatment + Type +
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

m3.T2.d15N <- gam(d15N ~
            s(plant.mass..g), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")


AIC.d15N.T2<-AIC(m1.T2.d15N, m2.T2.d15N, m3.T2.d15N)
## additive model best fit

summary(m1.T2.d15N)
anova.gam(m1.T2.d15N)
gam.check(m1.T2.d15N, rep=1000)
draw(m1.T2.d15N)
concrvity(m1.T2.d15N)
par(mfrow = c(1, 2))
plot(m1.T2.d15N, all.terms = TRUE, page=1)

# model predictions
d15N.diff.T2<-plot_difference(
  m1.T2.d15N,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
# model for smoothing
msmooth.T2.d15N<- gam(d15N ~ Treatment +
            s(plant.mass..g, by=Treatment), 
            subset = Time.point=="T2", data=topes.trt, method="REML", family="gaussian")

#plot for the model output on rawdata
d15N.T2.mod.plot<-
  plot_smooths(
  model = msmooth.T2.d15N,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=topes.trt[(topes.trt$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=d15N, color=Treatment, shape=Type)) +
  scale_shape_manual(name="Plankton", values = c(17, 16), 
                     labels = c(expression(paste("> 63"~mu,"m")), 
                                expression(paste("< 63"~mu,"m")))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 250)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

d15N.model<- plot_grid(d15N.T1.mod.plot + theme(legend.position = "none"), 
                    d15N.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend.mix, rel_widths = c(8,8,3), ncol=3)

ggsave("figures/d15N.model.pdf", encod="MacRoman", height=5, width=10)

################

mod.d15N<-rep(c("Treatment + Type + s(plant.mass..g, by=Treatment)", 
              "Treatment + Type + s(plant.mass..g)",
              "s(plant.mass..g)"), times=2)

mod.d15N.df<- data.frame(mod.d15N)
AIC.d15N<-bind_rows(AIC.d15N.T1, AIC.d15N.T2)

AIC.d15N.mod<-cbind(mod.d15N.df, AIC.d15N)
write.csv(AIC.d15N.mod, "output/AIC models/AIC.d15N.mod.csv")


#### figure for supplement
# legend
extract.legend.mix <- get_legend(
  # create some space to the left of the legend
  d15N.T2.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


## combine 
sage.mix.model.d15N<- plot_grid(d15N.T1.mod.plot + theme(legend.position = "none"), 
                    d15N.T2.mod.plot + theme(legend.position = "none"),
                    extract.legend.mix,
                    rel_widths = c(8,8,3), ncol=3, labels=c('A', 'B', ''), label_size=8)
sage.mix.model.d15N

ggsave("figures/Isotope.mixmodel.d15N.pdf", encod="MacRoman", height=4, width=8)

Greenhouse Gas Data

Samples for carbon dioxide (CO2) and methane (CH4) greenhouse gasses were collected from each tank on Days-0, 10, 31 and 59 of the experiment using the headspace method. Background concentrations of CO2 and CH4 in ambient air were collected at each sampling day by collecting 12 mL of air in evacuated Exetainers.

We use a script adapted from Koschorreck et al. 2021, accessed via github https://github.com/icra/headspace. This script produces a function that will run through processing of CO2 and CH4 ppm (raw data) and produce outputs in different units (ppm, molar, atm). For CO2, due to carbonate chemistry and the dissolution of CO2 in water, a simple headspace technique is not advised. Using the script of Koschorreck et al. 2021, we are able to account for CO2 dissolved in water and in other forms (bicarbonate and carbonate), however, an estimate or measurement of total alkalinity is required (see assumptions below).

####### ####### ####### ####### ####### ####### 
####### run function to generate outputs for the headspace technique
####### then load in the GHG data as raw ppm and run script 
####### ####### ####### ####### ####### ####### 


### Script adapted from Koschorreck et al. 2021, accessed via github https://github.com/icra/headspace

# Koschorreck M, Prairie YT, Kim J, Marcé R. 2021. Technical note: CO2 is not like CH4 – limits of and corrections to the headspace method to analyse pCO2 in fresh water. Biogeosciences  18:1619–1627. DOI: 10.5194/bg-18-1619-2021.

# script here will only calculate for freshwater system (option 1)
# for info on option 2 and 3 (brackish and marine waters), see Koschorreck et al. 2021

# below are annotations provided by the Dr. Marcé and colleagues

#####################################################################
# Rheadspace.R
#
# R function to calculate pCO2 in a water sample (micro-atm) using a complete headspace method accounting for the carbonate equilibrium in the equilibration vessel. 

# Authors: Rafael Marcé (Catalan Institute for Water Research - ICRA)
#          Jihyeon Kim (Université du Québec à Montréal - UQAM) 
#          Yves T. Prairie (Université du Québec à Montréal - UQAM)
#
# Contact information: Rafael Marcé (rmarce@icra.cat)
#####################################################################

# INPUT: 
#       You can either input a vector of 11 values for solving a single sample or a data frame of 11
#       columns and an arbitrary number of rows for batch processing of several samples.
#       
#       If supplying a vector for one sample, the vector should contain (in this order):
# 
#         1. The ID of the sample (arbitrary test, e.g."Sample_1")
#         2. mCO2 (ppmv) of the headspace "before" equilibration (e.g., zero for nitrogen)
#         3. mCO2 (ppmv) of the headspace "after" equilibration (e.g., as measured by a GC)
#         4. In situ (field) water temperature in degrees celsius
#         5. Water temperature after equilibration in degree celsius
#         6. Alkalinity (micro eq/L) of the water sample
#         7. Volume of gas in the headspace vessel (mL)
#         8. Volume of water in the headspace vessel (mL)
#         9. Barometric pressure at field conditions in kPa. 101.325 kPa = 1 atm 
#        10. Set of constants for carbonate equilibrium calculations 
                #(1=Freshwater, Millero 1979; 2=Estuarine, Millero 2010; 3=Marine, Dickson et al 2007) 
#        11. Salinity (PSU) # Set to zero if option in 10 is set to 1.
#
###########
#       If supplying a data frame, you can build it importing of a csv file
#       Example: dataset <- read.csv("R_test_data.csv")
#       The first row of this file must contain column names, then one row for each sample to be solved.
#       The columns names must be:
# 
#         1. Sample.ID 
#         2. HS.mCO2.before
#         3. HS.mCO2.after
#         4. Temp.insitu
#         5. Temp.equil
#         6. Alkalinity.measured
#         7. Volume.gas
#         8. Volume.water
#         9. Bar.pressure
#        10. Constants
#        11. Salinity
#
#       For the different samples, values must be as follows:
#
#         Sample.ID #User defined text
#         HS.mCO2.before #the pCO2 (ppmv) of the headspace "before" equilibration (e.g. zero for nitrogen)
#         HS.mCO2.after #the measured pCO2 (ppmv) of the headspace "after" equilibration
#         Temp.insitu #in situ (field) water temperature in degrees celsius
#         Temp.equil #the water temperature after equilibration in degree celsius
#         Alkalinity.measured #Total alkalinity (micro eq/L) of the water sample
#         Volume.gas #Volume of gas in the headspace vessel (mL)
#         Volume.water #Volume of water in the headspace vessel (mL)
#         Bar.pressure #Barometric pressure at field conditions in kPa. 101.325 kPa = 1 atm   
#         Constants #Set of constants for carbonate equilibrium calculations (1=Freshwater; 2=Estuarine; 3=Marine) 
#         Salinity # Salinity in PSU, set to zero if option in 10 is set to 1.
#
#
# EXAMPLE OF USE:
#  source("Rheadspace.R")
# 
#  pCO2 <- Rheadspace("Sample_1",0,80,20,25,1050,30,30,101.325,1,0)
# 
#  dataset <- read.csv("R_test_data.csv")
#  pCO2 <- Rheadspace(dataset)
#
# OUTPUT: a data frame containing:
#      1. Sample IDs
#      2. mCO2 complete headspace (ppmv) # mCO2 calculated using the complete headspace method accounting for the carbonate equilibrium
#      3. pCO2 complete headspace (micro-atm) # pCO2 calculated using the complete headspace method accounting for the carbonate equilibrium
#      4. CO2 concentration complete headspace (micro-mol/L) # CO2 concentration calculated using the complete headspace method accounting for the carbonate equilibrium
#      5. pH  # pH calculated for the sanple at in situ field conditions (using the complete headspace method)
#      6. mCO2 simple headspace (ppmv)  # mCO2 calculated using the simple headspace method NOT accounting for the carbonate equilibrium
#      7. pCO2 simple headspace (micro-atm)  # pCO2 calculated using the simple headspace method NOT accounting for the carbonate equilibrium
#      8. CO2 concentration simole headspace (micro-mol/L) # CO2 concentration calculated using the simple headspace method NOT accounting for the carbonate equilibrium
#      9. % error # error associated with using the simple headspace calculation
#
#
# REFERENCES
#
#  Dickson, A.G & J.P Riley (1979). The estimation of acid dissociation constants in sea-water media
#  from potentiometric titrations with strong base. II. The dissociation of phosphoric acid,
#  Marine Chemistry, 7(2), 101-109.  
#
#  Dickson, A. G., Sabine, C. L., and Christian, J. R. (2007). Guide to best practices for
#  ocean CO2 measurements, PICES Special Publication 3, 191 pp.
#
#  Millero, F. (1979). The thermodynamics of the carbonate system in seawater,
#  Geochimica et Cosmochimica Acta, 43(10), 1651-1661.  
#
#  Millero, F. (2010). Carbonate constants for estuarine waters,
#  Marine and Freshwater Research, 61(2), 139.
#
#  Orr, J. C., Epitalon, J.-M., and Gattuso, J.-P. (2015). Comparison of ten packages that compute
#  ocean carbonate chemistry, Biogeosciences, 12, 1483–1510. 
#
#  Weiss, R.F. (1974). Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
#  Marine Chemistry, 2, 203-215.
# 

#####################################################################
 
## THE FUNCTION 
## modified by C Wall for only freshwater system (salinity = 0, constants = 1)
## also modified to run for CO2 and CH4 from a common input csv dataframe

# first, the CO2 function

Rheadspace.CO2 <-  function(...){
  arguments <- list(...)
  
  # test arguments and initialize variables
  if (is.data.frame(arguments[[1]])) {
    input.table=arguments[[1]]
    if (dim(input.table)[2]!=11){
      stop("You should input a data frame with 11 columns. See the readme file or comments in the function", call.=FALSE)
    }else{
      Sample.ID = as.character(input.table$Sample.ID)
      mCO2_headspace = input.table$HS.mCO2.before #the mCO2 (ppmv) of the headspace "before" equilibration
      mCO2_eq = input.table$HS.mCO2.after #the measured mCO2 (ppmv) of the headspace "after" equilibration
      temp_insitu = input.table$Temp.insitu #in situ water temperature in degrees celsius
      temp_eq = input.table$Temp.equil #the water temperature after equilibration in degree celsius
      alk = input.table$Alkalinity.measured #Total alkalinity (micro eq/L) of the water sample
      vol_gas = input.table$Volume.gas #Volume of gas in the headspace vessel (mL)
      vol_water = input.table$Volume.water #Volume of water in the headspace vessel (mL)   
      Bar.pressure = input.table$Bar.pressure #Barometric pressure at field conditions in kPa. 101.325 kPa = 1 atm   
      c_constants = input.table$Constants #Constants for carbonate equilibrium (1=Freshwater; 2=Estuarine; 3=Marine) 
      Salinity = input.table$Salinity #Salinity in PSU. Set to zero if Constants = 1
      } 
  } else if (length(arguments)==11) {
    Sample.ID = as.character(arguments[[1]])
    mCO2_headspace = arguments[[2]] #the mCO2 (ppmv) of the headspace "before" equilibration
    mCO2_eq = arguments[[3]] #the measured mCO2 (ppmv) of the headspace "after" equilibration
    temp_insitu = arguments[[4]] #in situ water temperature in degrees celsius
    temp_eq = arguments[[5]] #the water temperature after equilibration in degree celsius
    alk = arguments[[6]] #Total alkalinity (micro eq/L) of the water sample
    vol_gas = arguments[[7]] #Volume of gas in the headspace vessel (mL)
    vol_water = arguments[[8]] #Volume of water in the headspace vessel (mL)   
    Bar.pressure = arguments[[9]] #Barometric pressure at field conditions in kPa. 101.325 kPa = 1 atm   
    c_constants = arguments[[10]] #Constants for carbonate equilibrium (1=Freshwater; 2=Estuarine; 3=Marine) 
    Salinity = arguments[[11]] #Salinity in PSU. Set to zero if Constants = 1
  } else {
    stop("You should input either a data frame or a vector of 11 values. See the readme file or comments in the function", call.=FALSE)
  }
  
  #initialization of variables -- this will be your output file of 9 columns
  pCO2_orig <- data.frame(matrix(NA,length(mCO2_headspace),11))
  
  names(pCO2_orig) <- c("Sample.ID","mCO2.complete.headspace..ppmv","pCO2.complete.headspace..micro.atm", "CO2 concentration.complete.headspace..micro.mol.L", "pH", "mCO2.simple.headspace..ppmv", "pCO2.simple.headspace..micro.atm.", "CO2.concentration.simple.headspace..micro.mol.L", "method.percent.error", "expected.Kh.of.CO2", "percent.excess.CO2")
 
  R <- 0.082057338 #L atm K-1 mol-1, for ideal gas law PV=nRT
  
  #the function uniroot cannot handle vectors, so we need a loop
  for (i in 1:length(mCO2_headspace)){ 
    
    AT = alk[i]*(1e-6) #conversion of TA in mEq/L to mol/L
    
    ##### Constants of the carbonate equilibrium
    # Kw = the dissociation constant of H2O into H+ and OH-
    # Kh = the solubility of CO2 in water - equilibration conditions
    # Kh2 = the solubility of CO2 in water - in situ field conditions
    # K1 = the equilibrium constant between CO2 and HCO3-
    # K2 = the equilibrium constant between HCO3- and CO3 2-
    
    # Solubility coefficients from Weiss (1974) with Sal=0 for freshwater option, using mols/L atm-1
    # Dissociation of water from Dickson and Riley (1979)
    
      #Millero, F. (1979). The thermodynamics of the carbonate system in seawater
      #Geochimica et Cosmochimica Acta 43(10), 1651 1661.  
      K1=10^-(-126.34048+6320.813/(temp_eq[i]+273.15)+19.568224*log(temp_eq[i]+273.15))
      K2=10^-(-90.18333+5143.692/(temp_eq[i]+273.15)+14.613358*log(temp_eq[i]+273.15))
      
      Kw = exp(148.9652-13847.26/(temp_eq[i]+273.15)-23.6521*log(273.15+temp_eq[i]))
      Kh = 10^((-58.0931+90.5069*(100/(273.15+temp_eq[i]))+22.294*log((273.15+temp_eq[i])/100))/log(10)) # mol/L/atm equilibration conditions
      Kh2 = 10^((-58.0931+90.5069*(100/(273.15+temp_insitu[i]))+22.294*log((273.15+temp_insitu[i])/100))/log(10)) # mol/L/atm original conditions
      
      
    HS.ratio <- vol_gas[i]/vol_water[i] #Headspace ratio (=vol of gas/vol of water)
    
    #The following calculations assume 1 atm, this is corrected later for measured pressure in the field
    #DIC at equilibrium
    co2 <- Kh * mCO2_eq[i]/1000000
    h_all <- polyroot(c(-(2*K1*K2*co2),-(co2*K1+Kw),AT,1))
    real<-Re(h_all)
    h <-real[which(real>0)]
    DIC_eq <- co2 * (1 + K1/h + K1 * K2/(h * h))
    
    #DIC in the original sample
    DIC_ori <- DIC_eq + (mCO2_eq[i] - mCO2_headspace[i])/1000000/(R*(temp_eq[i]+273.15))*HS.ratio
    
    #pCO2 in the original sample
    h_all <- polyroot(c(-(K1*K2*Kw),K1*K2*AT-K1*Kw-2*DIC_ori*K1*K2,AT*K1-Kw+K1*K2-DIC_ori*K1,AT+K1,1))
    real<-Re(h_all)
    h <-real[which(real>0)]
    
    co2 <- h* (DIC_ori * h * K1/(h * h + K1 * h + K1 * K2)) / K1
    
    pCO2_orig[i,1] <- as.character(Sample.ID[i])
    pCO2_orig[i,2] <- co2/Kh2*1000000
    pCO2_orig[i,3] <- pCO2_orig[i,2]*Bar.pressure[i]/101.325
    pCO2_orig[i,4] <- co2*1000000
    pCO2_orig[i,5] <- -log10( h )
    
    
    # Calculation not accounting for alkalinity effects and associated error, this is the "simple" method
    
    #concentration and total mass in the water sample assuming ideal gas from the pCO2 measured at the headspace
    CO2_solution <- mCO2_eq[i]/1000000*Kh #mol/L
    CO2_solution_mass <- CO2_solution * vol_water[i]/1000 #mol
    
    #mass of CO2 in the measured headspace
    final_C_headspace_mass <- mCO2_eq[i]/1000000*(vol_gas[i]/1000) / (R * (temp_eq[i]+273.15)) #mol
    
    mols_headspace <- mCO2_headspace[i]/1000000*(vol_gas[i]/1000)/(R * (temp_eq[i]+273.15)) #mol PV / RT = n
    
    #implication: mass, concentration, and partial pressure of CO2 in the original sample (amount in sample and headspace after equilibration minus original mass in the headspace)
    Sample_CO2_mass <- CO2_solution_mass + final_C_headspace_mass - mols_headspace #mol
    Sample_CO2_conc <- Sample_CO2_mass/(vol_water[i]/1000) #mol/L
    pCO2_orig[i,6] <- Sample_CO2_conc/Kh2*1000000 #ppmv
    pCO2_orig[i,7] <- pCO2_orig[i,6]*Bar.pressure[i]/101.325 # micro-atm
    pCO2_orig[i,8] <- Sample_CO2_conc*1000000 # micro-mol/L
    #calculation of the error
    pCO2_orig[i,9] <- (pCO2_orig[i,6]-pCO2_orig[i,2])/pCO2_orig[i,2] *100  #%
    
    # calculate the percent excess
    expected.uM<- Kh*(mCO2_headspace[i]) # kH in mol/L * ppm background 
    # note: conversion above to mol would be mCO2_headspace[i]/10^6, but then * by 10^6 to umol, so 10^6 cancels.
    pCO2_orig[i,10]<- expected.uM #umol/L
    measured.uM<-co2*1000000  # where this value is the pCO2_orig[i,4] at umol/L
    pCO2_orig[i,11]<- (measured.uM / expected.uM) *100 # percent excess
  }
  
  
  return(pCO2_orig) #Output data frame
  
}



##########################################################################################
## Methane function
##########################################################################################
Rheadspace.CH4 <-  function(...){
  arguments <- list(...)
  
  # test arguments and initialize variables
  if (is.data.frame(arguments[[1]])) {
    input.table=arguments[[1]]
    if (dim(input.table)[2]!=11){
      stop("You should input a data frame with 11 columns. See the readme file or comments in the function", call.=FALSE)
    }else{
      Sample.ID = as.character(input.table$Sample.ID)
      mCH4_headspace = input.table$HS.mCH4.before #the mCH4 (ppmv) of the headspace "before" equilibration
      mCH4_eq = input.table$HS.mCH4.after #the measured mCH4 (ppmv) of the headspace "after" equilibration
      temp_insitu = input.table$Temp.insitu #in situ water temperature in degrees celsius
      temp_eq = input.table$Temp.equil #the water temperature after equilibration in degree celsius
      alk = input.table$Alkalinity.measured #Total alkalinity (micro eq/L) of the water sample
      vol_gas = input.table$Volume.gas #Volume of gas in the headspace vessel (mL)
      vol_water = input.table$Volume.water #Volume of water in the headspace vessel (mL)   
      Bar.pressure = input.table$Bar.pressure #Barometric pressure at field conditions in kPa. 101.325 kPa = 1 atm   
      c_constants = input.table$Constants #Constants for carbonate equilibrium (1=Freshwater; 2=Estuarine; 3=Marine) 
      Salinity = input.table$Salinity #Salinity in PSU. Set to zero if Constants = 1
    } 
  } else if (length(arguments)==11) {
    Sample.ID = as.character(arguments[[1]])
    mCH4_headspace = arguments[[2]] #the mCH4 (ppmv) of the headspace "before" equilibration
    mCH4_eq = arguments[[3]] #the measured mCH4 (ppmv) of the headspace "after" equilibration
    temp_insitu = arguments[[4]] #in situ water temperature in degrees celsius
    temp_eq = arguments[[5]] #the water temperature after equilibration in degree celsius
    alk = arguments[[6]] #Total alkalinity (micro eq/L) of the water sample
    vol_gas = arguments[[7]] #Volume of gas in the headspace vessel (mL)
    vol_water = arguments[[8]] #Volume of water in the headspace vessel (mL)   
    Bar.pressure = arguments[[9]] #Barometric pressure at field conditions in kPa. 101.325 kPa = 1 atm   
    c_constants = arguments[[10]] #Constants for carbonate equilibrium (1=Freshwater; 2=Estuarine; 3=Marine) 
    Salinity = arguments[[11]] #Salinity in PSU. Set to zero if Constants = 1
  } else {
    stop("You should input either a data frame or a vector of 11 values. See the readme file or comments in the function", call.=FALSE)
  }
  
  
  #initialization of variables
  pCH4_orig <- data.frame(matrix(NA,length(mCH4_headspace),6))
  names(pCH4_orig) <- c("Sample.ID", "mCH4.simple.headspace..ppmv", "pCH4.simple.headspace..micro.atm", 
                        "CH4.concentration.simple.headspace..nano.mol.L", "expected.Kh.of.CH4", "percent.excess.CH4")
  
  R <- 0.082057338 #L atm K-1 mol-1
  
  #the function uniroot cannot handle vectors, so we need a loop
  for (i in 1:length(mCH4_headspace)){ 
    
    
    ## CH4 Solubility coefficients from Yamamoto et al (1976) with Sal=0 for freshwater option
    #Yamamoto S, Alcauskas JB, Crozier TE. 1976. Solubility of methane in distilled water and seawater. Journal of chemical and engineering data 21:78–80. DOI: 10.1021/je60068a029.
    
    Kh = exp(-67.1962+99.1624*(100/(273.15+temp_eq[i]))+27.9015*log((273.15+temp_eq[i])/100))
    # LCH4/LH2O/atm equilibration conditions
    
    Kh2 = exp(-67.1962+99.1624*(100/(273.15+temp_insitu[i]))+27.9015*log((273.15+temp_insitu[i])/100))
    # this isn't applied fully, but if post-equilibration temp different, or measured, this would be useful
    
    #converting coefficient from L/L to mol CH4 /L H2O using solubility/RT, equilibration conditions
    Kh.moles.L<-Kh/(0.082057338*(273.15+temp_eq[i]))
    
    Kh2.moles.L<-Kh2/(0.082057338*(273.15+temp_insitu[i])) #original temp
    
    
    #Calculation headspace 
    
    #concentration and total mass in the water sample assuming ideal gas from the pCH4 measured at the headspace
    CH4_solution <- mCH4_eq[i]/1000000*Kh.moles.L # ppm headspace mol/L
    CH4_solution_mass <- CH4_solution * (vol_water[i]/1000) # mass in mol
    
    #mass of CH4 in the measured headspace
    final_CH4_headspace_mass <- mCH4_eq[i]/1000000*(vol_gas[i]/1000) / (R * (temp_eq[i]+273.15)) #mol
    
    molsCH4_headspace <- mCH4_headspace[i]/1000000*(vol_gas[i]/1000) / (R * (temp_eq[i]+273.15)) #mol PV / RT = n
    
    #implication: mass, concentration, and partial pressure of CH4 in the original sample (amount in sample and headspace after equilibration minus original mass in the headspace)
    
    Sample_CH4_mass <- CH4_solution_mass + final_CH4_headspace_mass - molsCH4_headspace #mol
    Sample_CH4_conc <- Sample_CH4_mass/(vol_water[i]/1000) #mol/L
    CH4.sample..ppm <- Sample_CH4_conc/Kh2.moles.L*1000000 #ppmv
    CH4.sample..uatm <- CH4.sample..ppm*(Bar.pressure[i]/101.325) # micro-atm
    CH4..nmol.L <- Sample_CH4_conc*10^9 # nano-mol/L
    
    # build output
    pCH4_orig[i,1] <- as.character(Sample.ID[i])
    pCH4_orig[i,2] <- CH4.sample..ppm
    pCH4_orig[i,3] <- CH4.sample..uatm
    pCH4_orig[i,4] <- CH4..nmol.L
    
    # calculate the percent excess
    expected.uM<- Kh.moles.L*(mCH4_headspace[i]) # kH in mol/L * ppm background 
    # note: conversion above to mol would be mCO2_headspace[i]/10^6, but then * by 10^6 to umol, so 10^6 cancels.
    expected.nM<-expected.uM*1000 #expected GHG as nmol/L, converted to nmol from above
    pCH4_orig[i,5]<- expected.nM
    measured.nM<- CH4..nmol.L  # where this value is the pCH4_orig[i,4] at nmol/L
    pCH4_orig[i,6]<- (measured.nM / expected.nM) *100 # percent excess > ambient
  }
  
  return(pCH4_orig) #Output data frame
  
}

Important assumptions in our data for GHG calculations.

  1. we estimate an average alkalinity for tanks based on Miramar tap/drinking water. Data from City of San Diego shows Miramar water with average of = 112 ppm (equivalent to mgCaCO3/L). Code requires units to be in mEq/L, so we convert to Eq/L by dividing mgCaCO3 by 50 (112/50= 2.24, to mEq/L = 2240).

  2. While salinity was determined here, these values are <1 unit salinty (g/L or ppt) and were only estimated for all tanks using a mean temperature on the day of collection. Therefore, we run the script without consideration for salinity, and this is reflected in our solubility equations from Weis et al. (1974) and Yamaoto et al. (1976), for CO2 and CH4, respectively.

  3. In the original version of this script (Koschorreck et al. 2021, accessed via github https://github.com/icra/headspace ) there are 3 options (ie., constants) for freshwater (=1), brackish/estuary (=2), marine (=3). We only use option 1 and have removed code for the function. See the original code for further info.

  4. We have added a section of script to calculate a percent excess for both GHGs. Here, at quilibrium, the expected concentration of GHGs in water is equal to ambient air (or “before shaking”) in the syringe. The observed or measured GHG is what is the concentration of gases in the water based on the headspace technique and script outputs. The ratio of these (measured / expected * 100) is the percent excess or saturation of GHGs relative to the atmosphere.

References

  • Koschorreck M, Prairie YT, Kim J, Marcé R. 2021. Technical note: CO2 is not like CH4 – limits of and corrections to the headspace method to analyse pCO2 in fresh water. Biogeosciences 18:1619–1627. DOI: 10.5194/bg-18-1619-2021.

  • Yamamoto S, Alcauskas JB, Crozier TE. 1976. Solubility of methane in distilled water and seawater. Journal of chemical and engineering data 21:78–80. DOI: 10.1021/je60068a029.

  • Weiss RF. 1974. Carbon dioxide in water and seawater: the solubility of a non-ideal gas. Marine chemistry 2:203–215. DOI: 10.1016/0304-4203(74)90015-2.

##########################################################################################
# load the data
#library("dplyr")
dataset <- read.csv("data/GH.gases/Pyro_CO2_CH4_input.csv")

# subset the metadata here to yield 11 columns that will be used to run CO2 and CH4 analysis
metadata.cols<-dataset %>%
  select(Time.point, Date, Treatment, plant.mass..g, Tank)

# subset the CO2 datset
dataset.CO2<-dataset %>%
  select(Sample.ID, HS.mCO2.before, HS.mCO2.after, Temp.insitu, Temp.equil, Alkalinity.measured,
         Volume.gas, Volume.water, Bar.pressure, Constants, Salinity)

# subset the CH4 datset
dataset.CH4<-dataset %>%
  select(Sample.ID, HS.mCH4.before, HS.mCH4.after, Temp.insitu, Temp.equil, Alkalinity.measured,
         Volume.gas, Volume.water, Bar.pressure, Constants, Salinity)

##########################################################################################

# run the function
pCO2 <- Rheadspace.CO2(dataset.CO2)
pCH4 <- Rheadspace.CH4(dataset.CH4)

# combine output with metadata to generate 2 dataframes
pCO2.output<-as.data.frame(c(metadata.cols, pCO2))
pCH4.output<-as.data.frame(c(metadata.cols, pCH4))

# export the data
#write.csv(pCO2.output, "data/GH.gases/pCO2_calc.output.csv")
#write.csv(pCH4.output, "data/GH.gases/pCH4_calc.output.csv")
########################################################################################## 

# combine to be dataframe for analysis

GHG.calc<-as.data.frame(cbind(metadata.cols, 
                pCO2.output$mCO2.complete.headspace..ppmv, 
                pCO2.output$CO2.concentration.complete.headspace..micro.mol.L, 
                pCO2.output$expected.Kh.of.CO2, 
                pCO2.output$percent.excess.CO2,
                
                pCH4.output$mCH4.simple.headspace..ppmv,
                pCH4.output$CH4.concentration.simple.headspace..nano.mol.L,
                pCH4.output$expected.Kh.of.CH4, 
                pCH4.output$percent.excess.CH4))

# rename columns as  new = old
GHG<- GHG.calc %>%
  rename(CO2.ppm =  "pCO2.output$mCO2.complete.headspace..ppmv",
         CO2.uM = "pCO2.output$CO2.concentration.complete.headspace..micro.mol.L",
         CO2.expect.uM = "pCO2.output$expected.Kh.of.CO2",
         perc.excess.CO2 = "pCO2.output$percent.excess.CO2",
         
         CH4.ppm = "pCH4.output$mCH4.simple.headspace..ppmv",
         CH4.nM = "pCH4.output$CH4.concentration.simple.headspace..nano.mol.L",
         CH4.expect.nM = "pCH4.output$expected.Kh.of.CH4",
         perc.excess.CH4 = "pCH4.output$percent.excess.CH4")

write.csv(GHG, "data/GH.gases/GHG.calc.output.csv")             

#####################
### we will be using the 'GHG' dataframe from above

# if reloading the data...
#GHG<-read.csv("data/GH.gases/GHG.calc.output.csv", header=TRUE)
#GHG<-na.omit(GHG) # remove NAs, this removes the blanks
#GHG<-GHG[, -1] # removes the "X column" which comes in when you load in the exported data above


# set structure
make.fac<-c("Time.point", "Treatment", "Tank")
GHG[make.fac] <- lapply(GHG[make.fac], factor) # make all these factors
GHG$plant.mass..g<-as.numeric(GHG$plant.mass..g)

# make a CO2 and CH4 dataframe
CO2.GHG<- GHG %>% 
  select(Time.point, Date, Treatment, plant.mass..g, Tank, CO2.ppm, CO2.uM, CO2.expect.uM, perc.excess.CO2)

CH4.GHG<- GHG %>% 
  select(Time.point, Date, Treatment, plant.mass..g, Tank, CH4.ppm, CH4.nM, CH4.expect.nM, perc.excess.CH4)

# remove any values less than 0, a few early measures of CH4 were lower than atmosphere
CH4.GHG<-CH4.GHG[!(CH4.GHG$CH4.nM < 0 ),] 

Now read in the output GHG data and run through model fitting and generate plots, statistical tests.


#### Make plots and run models for CO2

CO2<-na.omit(CO2.GHG) # dataframe in ppm and uM of CO2

CO2.expect<-aggregate(CO2.expect.uM~ Time.point, data=CO2, FUN=mean) # this is for the mean of the expected CO2

############# GHG plots and modesl
#### CO2
#######-- T0

m1.T0.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T0", data = CO2,
            method="REML", family="gaussian")

m2.T0.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T0", data = CO2,
            method="REML", family="gaussian")

m3.T0.CO2<- gam(CO2.uM ~
            s(plant.mass..g), subset = Time.point=="T0", data = CO2,
            method="REML", family="gaussian")

CO2.T0.AIC<-AIC(m1.T0.CO2,m2.T0.CO2, m3.T0.CO2)
#factor smooth best

summary(m1.T0.CO2)
anova.gam(m1.T0.CO2)
gam.check(m1.T0.CO2, rep=1000)
draw(m1.T0.CO2)
concrvity(m1.T0.CO2)
par(mfrow = c(1, 2))
plot(m1.T0.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T0<-plot_difference(
  m1.T0.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T0.mod.plot<-
  plot_smooths(
  model = m1.T0.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=CO2.uM, color=Treatment)) +
  geom_hline(yintercept= CO2.expect[1,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-0") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T0.mod.plot

#### CO2 
#######-- T1

m1.T1.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T1", data = CO2,
            method="REML", family="gaussian")

m2.T1.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T1", data = CO2,
            method="REML", family="gaussian")

m3.T1.CO2<- gam(CO2.uM ~
            s(plant.mass..g), subset = Time.point=="T1", data = CO2,
            method="REML", family="gaussian")


CO2.T1.AIC<-AIC(m1.T1.CO2,m2.T1.CO2, m3.T1.CO2)
#factor smooth best

summary(m1.T1.CO2)
anova.gam(m1.T1.CO2)
gam.check(m1.T1.CO2, rep=1000)
draw(m1.T1.CO2)
concrvity(m1.T1.CO2)
par(mfrow = c(1, 2))
plot(m1.T1.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T1<-plot_difference(
  m1.T1.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T1.mod.plot<-
  plot_smooths(
  model = m1.T1.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=CO2.uM, color=Treatment)) +
  geom_hline(yintercept= CO2.expect[2,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T1.mod.plot

#### CO2
#######-- T2

m1.T2.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T2", data = CO2,
            method="REML", family="gaussian")

m2.T2.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T2", data = CO2,
            method="REML", family="gaussian")

m3.T2.CO2<- gam(CO2.uM ~
            s(plant.mass..g), subset = Time.point=="T2", data = CO2,
            method="REML", family="gaussian")

CO2.T2.AIC<-AIC(m1.T2.CO2, m2.T2.CO2, m3.T2.CO2)
# smooth by factor best

summary(m1.T2.CO2)
anova.gam(m1.T2.CO2)
gam.check(m1.T2.CO2, rep=1000)
draw(m1.T2.CO2)
concrvity(m1.T2.CO2)
par(mfrow = c(1, 2))
plot(m1.T2.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T2<-plot_difference(
  m1.T2.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T2.mod.plot<-
  plot_smooths(
  model = m1.T2.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=CO2.uM, color=Treatment)) +
  geom_hline(yintercept= CO2.expect[3,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T2.mod.plot


#### CO2 
#######-- T3

m1.T3.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T3", data = CO2,
            method="REML", family="gaussian")

m2.T3.CO2<- gam(CO2.uM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T3", data = CO2,
            method="REML", family="gaussian")

m3.T3.CO2<- gam(CO2.uM ~
            s(plant.mass..g), subset = Time.point=="T3", data = CO2,
            method="REML", family="gaussian")

CO2.T3.AIC<-AIC(m1.T3.CO2, m2.T3.CO2, m3.T3.CO2)
# smooth by factor best

summary(m1.T3.CO2)
anova.gam(m1.T3.CO2)
gam.check(m1.T3.CO2, rep=1000)
draw(m1.T3.CO2)
concrvity(m1.T3.CO2)
par(mfrow = c(1, 2))
plot(m1.T3.CO2, all.terms = TRUE, page=1)

# model predictions
CO2.diff.T3<-plot_difference(
  m1.T3.CO2,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CO2.T3.mod.plot<-
  plot_smooths(
  model = m1.T3.CO2,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CO2[(CO2$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=CO2.uM, color=Treatment)) +
  geom_hline(yintercept= CO2.expect[4,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CO"[2], ~(mu*M), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-3") +
  coord_cartesian(ylim=c(0, 400)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CO2.T3.mod.plot

# models and raw data
CO2.model.plots<-plot_grid(
  CO2.T0.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-0"),
  CO2.T1.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-10"),
  CO2.T2.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-31"),
  CO2.T3.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-59"),
  extract.legend, 
  rel_widths = c(8,8,8,8,3), ncol=5)

### model differences
CO2.plot.diff<-plot_grid(
  CO2.diff.T0+ ggtitle("CO2-T0"),
  CO2.diff.T1+ ggtitle("Day-10"),
  CO2.diff.T2+ ggtitle("Day-31"),
  CO2.diff.T3+ ggtitle("Day-59"),
  rel_widths = c(8,8,8,8), ncol=4)

## AIC table
mod.ghg<-rep(c("Treatment +s(plant.mass..g, by=Treatment)", 
              "Treatment + s(plant.mass..g)", 
              "s(plant.mass..g)"), times=4)

mod.ghg.df<- data.frame(mod.ghg)

AIC.CO2<-bind_rows(CO2.T0.AIC, CO2.T1.AIC, CO2.T2.AIC, CO2.T3.AIC)
AIC.CO2.mod<-cbind(mod.ghg.df, AIC.CO2)
write.csv(AIC.CO2.mod, "output/AIC models/AIC.CO2.mod.csv")

Table: Results for CO2 across burned/unburned treatments at Day-0.

anova.gam(m1.T0.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CO2.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F  p-value
## Treatment  1 14.62 0.000869
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value
## s(plant.mass..g):Treatmentburned   1.000  1.000 11.505 0.00251
## s(plant.mass..g):Treatmentunburned 3.966  4.815  3.342 0.02586

Table: Results for CO2 across burned/unburned treatments at Day-10.

anova.gam(m1.T1.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CO2.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df    F p-value
## Treatment  1 9.68 0.00494
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.070  2.556 166.0  <2e-16
## s(plant.mass..g):Treatmentunburned 3.062  3.748 150.6  <2e-16

Table: Results for CO2 across burned/unburned treatments at Day-31.

anova.gam(m1.T2.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CO2.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.344   0.563
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   2.661  3.267 51.80  <2e-16
## s(plant.mass..g):Treatmentunburned 3.495  4.261 23.07  <2e-16

Table: Results for CO2 across burned/unburned treatments at Day-59.

anova.gam(m1.T3.CO2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CO2.uM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 2.476   0.129
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value
## s(plant.mass..g):Treatmentburned   2.808  3.443 13.27 1.70e-05
## s(plant.mass..g):Treatmentunburned 1.000  1.000 29.70 1.45e-05

Now make plots and run models for methane.

CH4<- na.omit(CH4.GHG) # CH4 data in ppm and nM

CH4.expect<-aggregate(CH4.expect.nM~ Time.point, data=CH4, FUN=mean) # mean of expected CH4

m1.T0.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T0", data = CH4,
            method="REML", family="gaussian")

m2.T0.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T0", data = CH4,
            method="REML", family="gaussian")

m3.T0.CH4<- gam(CH4.nM ~
            s(plant.mass..g), subset = Time.point=="T0", data = CH4,
            method="REML", family="gaussian")

CH4.T0.AIC<-AIC(m1.T0.CH4, m2.T0.CH4, m3.T0.CH4)
#factor single smooth best

summary(m2.T0.CH4)
anova.gam(m2.T0.CH4)
gam.check(m2.T0.CH4, rep=1000)
draw(m2.T0.CH4)
concrvity(m2.T0.CH4)
par(mfrow = c(1, 2))
plot(m2.T0.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T0<-plot_difference(
  m2.T0.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T0.mod.plot<-
  plot_smooths(
  model = m2.T0.CH4,
  series = plant.mass..g,
  comparison=Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=CH4.nM, color=Treatment)) +
  geom_hline(yintercept= CH4.expect[1,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  ylab(expression(paste("CH"[4], ~(nM), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-0") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CH4.T0.mod.plot

#### CH4 
#######-- T1

m1.T1.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T1", data = CH4,
            method="REML", family="gaussian")

m2.T1.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T1", data = CH4,
            method="REML", family="gaussian")

m3.T1.CH4<- gam(CH4.nM ~
            s(plant.mass..g), subset = Time.point=="T1", data = CH4,
            method="REML", family="gaussian")

CH4.T1.AIC<-AIC(m1.T1.CH4, m2.T1.CH4, m3.T1.CH4)
#factor by smooth best

summary(m1.T1.CH4)
anova.gam(m1.T1.CH4)
gam.check(m1.T1.CH4, rep=1000)
draw(m1.T1.CH4)
concrvity(m1.T1.CH4)
par(mfrow = c(1, 2))
plot(m1.T1.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T1<-plot_difference(
  m1.T1.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T1.mod.plot<-
  plot_smooths(
  model = m1.T1.CH4,
  series = plant.mass..g,
  comparison = Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=CH4.nM, color=Treatment)) +
  geom_hline(yintercept= CH4.expect[2,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CH"[4], ~(nM), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-1") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CH4.T1.mod.plot

#### CH4
#######-- T2

m1.T2.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T2", data = CH4,
            method="REML", family="gaussian")

m2.T2.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T2", data = CH4,
            method="REML", family="gaussian")

m3.T2.CH4<- gam(CH4.nM ~
            s(plant.mass..g), subset = Time.point=="T2", data = CH4,
            method="REML", family="gaussian")

CH4.T2.AIC<-AIC(m1.T2.CH4, m2.T2.CH4, m3.T2.CH4)
# global best

summary(m3.T2.CH4)
anova.gam(m3.T2.CH4)
gam.check(m3.T2.CH4, rep=1000)
draw(m3.T2.CH4)
concrvity(m3.T2.CH4)
par(mfrow = c(1, 2))
plot(m3.T2.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T2<-plot_difference(
  m1.T2.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T2.mod.plot<-
  plot_smooths(
  model = m3.T2.CH4,
  series = plant.mass..g,
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=CH4.nM, color=Treatment)) +
  geom_hline(yintercept= CH4.expect[3,2], linetype="longdash", color = "gray") +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ylab(expression(paste("CH"[4], ~(nM), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-2") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))

CH4.T2.mod.plot


#### CH4 
#######-- T3

m1.T3.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g, by=Treatment), subset = Time.point=="T3", data = CH4,
            method="REML", family="gaussian")

m2.T3.CH4<- gam(CH4.nM ~ Treatment +
            s(plant.mass..g), subset = Time.point=="T3", data = CH4,
            method="REML", family="gaussian")

m3.T3.CH4<- gam(CH4.nM ~
            s(plant.mass..g), subset = Time.point=="T3", data = CH4,
            method="REML", family="gaussian")

CH4.T3.AIC<-AIC(m1.T3.CH4, m2.T3.CH4, m3.T3.CH4)
# global with treatment best

summary(m2.T3.CH4)
anova.gam(m2.T3.CH4)
gam.check(m2.T3.CH4, rep=1000)
draw(m2.T3.CH4)
concrvity(m2.T3.CH4)
par(mfrow = c(1, 2))
plot(m2.T3.CH4, all.terms = TRUE, page=1)

# model predictions
CH4.diff.T3<-plot_difference(
  m1.T3.CH4,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
#plot for the model output on rawdata
CH4.T3.mod.plot<-
  plot_smooths(
  model = m2.T3.CH4,
  series = plant.mass..g,
  comparison=Treatment
)  + theme(legend.position = "none") +
  geom_point(data=CH4[(CH4$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=CH4.nM, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= CH4.expect[4,2], linetype="longdash", color = "gray") +
  geom_line(aes(fill=Treatment, linetype=Treatment)) +
  ylab(expression(paste("CH"[4], ~(nM), sep=""))) +
  xlab("plant material (g)") +
  ggtitle("Time-3") +
  coord_cartesian(ylim=c(0, 50)) +
  Fig.formatting +
  theme(legend.key.size = unit(1,"line"))


CH4.model.plots<-plot_grid(
  CH4.T0.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-0"),
  CH4.T1.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-10"),
  CH4.T2.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-31"),
  CH4.T3.mod.plot+ theme(legend.position= "none") + ggtitle( "Day-59"), 
  extract.legend, 
  rel_widths = c(8,8,8,8,3), ncol=5)

### model differences
CH4.plot.diff<-plot_grid(
  CH4.diff.T0+ ggtitle("CH4-Day-0"),
  CH4.diff.T1+ ggtitle("Day-10"),
  CH4.diff.T2+ ggtitle("Day-31"),
  CH4.diff.T3+ ggtitle("Day-59"),
  rel_widths = c(8,8,8,8), ncol=4)


##### AIC table
AIC.CH4<-bind_rows(CH4.T0.AIC, CH4.T1.AIC, CH4.T2.AIC, CH4.T3.AIC)

AIC.CH4.mod<-cbind(mod.ghg.df, AIC.CH4)
write.csv(AIC.CH4.mod, "output/AIC models/AIC.CH4.mod.csv")

Table: Results for CH4 across burned/unburned treatments at Day-0.

anova.gam(m2.T0.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CH4.nM ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 2.038   0.166
## 
## Approximate significance of smooth terms:
##                  edf Ref.df     F p-value
## s(plant.mass..g)   1      1 0.718   0.405

Table: Results for CH4 across burned/unburned treatments at Day-10.

anova.gam(m1.T1.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CH4.nM ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.266   0.611
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.812  2.243 0.889 0.42790
## s(plant.mass..g):Treatmentunburned 3.346  4.086 6.531 0.00119

Table: Results for CH4 across burned/unburned treatments at Day-31.

anova.gam(m3.T2.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CH4.nM ~ s(plant.mass..g)
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    F p-value
## s(plant.mass..g) 1.000  1.001 0.19   0.667

Table: Results for CH4 across burned/unburned treatments at Day-59.

anova.gam(m2.T3.CH4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## CH4.nM ~ Treatment + s(plant.mass..g)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 3.646  0.0679
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F p-value
## s(plant.mass..g) 3.381  4.127 2.037   0.113

Compile the greenhouse gas plots for CO2 and CH4. Fit the data with model fitting and plot the model difference.

GHG.plots<-plot_grid(CO2.model.plots, CH4.model.plots, ncol=1, labels=c('A', 'B'), label_size=8)
GHG.plots
**Figure** Greenhouse gas concentration for (**A**) carbon dioxide (CO~2~) and (**B**) methane (CH~4~) at the beginning of the study before plant material was added (Day-0) and across three experimental time points. Horizontal gray lines indicate the atmospheric gas concentrations at each sampling day.  Lines represent best-fit generalized additive models (GAMs) with 95% confidence intervals.  Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

Figure Greenhouse gas concentration for (A) carbon dioxide (CO2) and (B) methane (CH4) at the beginning of the study before plant material was added (Day-0) and across three experimental time points. Horizontal gray lines indicate the atmospheric gas concentrations at each sampling day. Lines represent best-fit generalized additive models (GAMs) with 95% confidence intervals. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

ggsave("figures/GHG.molar.mod.plots.pdf", height=7, width=12)
GHG.plot.diff<-plot_grid(CO2.plot.diff, CH4.plot.diff, ncol=1, labels=c('A', 'B'), label_size=8)
GHG.plot.diff
**Figure** Model effects from GAMs with differences between smoothers for greenhouse gasses (**A**) carbon dioxide (CO~2~) and (**B**) methane (CH~4~) in tanks receiving burned and unburned plant material at the beginning of the experiment and during three experimental time points. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Figure Model effects from GAMs with differences between smoothers for greenhouse gasses (A) carbon dioxide (CO2) and (B) methane (CH4) in tanks receiving burned and unburned plant material at the beginning of the experiment and during three experimental time points. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

ggsave("figures/GHG.molar.mod.diff.pdf", height=7, width=12)

Percent excess GHG

Using the GHG data, now produce a "percent excess" to show how much more CO2 or CH4 (molar scale) is present in the mesocosms versus the atmosphere (= expected GHG concentration at equilibrium). This data is essentially the same shape and relationship as the model-fit data above, but with a more simplistic fitting via geom_smooth and a simple gam smoother function.

######### percent excess CO2

T0.CO2<-ggplot(data=CO2[(CO2$Time.point=="T0"),], aes(y=perc.excess.CO2, x=plant.mass..g, fill=Treatment))+
  stat_smooth(aes(y=perc.excess.CO2,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam") +
  ylab(expression(paste("% excess CO"[2], sep=""))) +
  coord_cartesian(ylim=c(0, 2000)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  ggtitle("Day-0")+
  theme_classic()+
  theme(legend.position= "none")

T1.CO2<-ggplot(data=CO2[(CO2$Time.point=="T1"),], aes(y=perc.excess.CO2, x=plant.mass..g, fill=Treatment))+
  stat_smooth(aes(y=perc.excess.CO2,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam", alpha=0.2) +
  ylab(expression(paste("% excess CO"[2], sep=""))) +
  coord_cartesian(ylim=c(0, 2000)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  ggtitle("Day-10")+
  theme_classic() +
  theme(legend.position= "none")

T2.CO2<-ggplot(data=CO2[(CO2$Time.point=="T2"),], aes(y=perc.excess.CO2, x=plant.mass..g, fill=Treatment))+
 stat_smooth(aes(y=perc.excess.CO2,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam", alpha=0.2) +
  ylab(expression(paste("% excess CO"[2], sep=""))) +
  coord_cartesian(ylim=c(0, 2000)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  ggtitle("Day-31")+
  theme_classic()+
  theme(legend.position= "none")

T3.CO2<-ggplot(data=CO2[(CO2$Time.point=="T3"),], aes(y=perc.excess.CO2, x=plant.mass..g, fill=Treatment))+
  stat_smooth(aes(y=perc.excess.CO2,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam", alpha=0.2) +
  ylab(expression(paste("% excess CO"[2], sep=""))) +
  coord_cartesian(ylim=c(0, 2000)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  ggtitle("Day-59")+
  theme_classic()+
  theme(legend.position= "none")


################## methane excess

T0.CH4<-ggplot(data=CH4[(CH4$Time.point=="T0"),], aes(y=perc.excess.CH4, x=plant.mass..g, fill=Treatment))+
  stat_smooth(aes(y=perc.excess.CH4,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam", alpha=0.2) +
  ylab(expression(paste("% excess CH"[4], sep=""))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  coord_cartesian(ylim=c(0, 1000)) +
  ggtitle("Day-0")+
  theme_classic()+
  theme(legend.position= "none")

T1.CH4<-ggplot(data=CH4[(CH4$Time.point=="T1"),], aes(y=perc.excess.CH4, x=plant.mass..g, fill=Treatment))+
 stat_smooth(aes(y=perc.excess.CH4,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam", alpha=0.2) +
  ylab(expression(paste("% excess CH"[4], sep=""))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  coord_cartesian(ylim=c(0, 1000)) +
  ggtitle("Day-10")+
  theme_classic() +
  theme(legend.position= "none")

T2.CH4<-ggplot(data=CH4[(CH4$Time.point=="T2"),], aes(y=perc.excess.CH4, x=plant.mass..g, fill=Treatment))+
  stat_smooth(aes(y=perc.excess.CH4,x=plant.mass..g, color=Treatment, linetype=Treatment),method="gam",alpha=0.2) +
  ylab(expression(paste("% excess CH"[4], sep=""))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  coord_cartesian(ylim=c(0, 1000)) +
  ggtitle("Day-31")+
  theme_classic()+
  theme(legend.position= "none")

T3.CH4<-ggplot(data=CH4[(CH4$Time.point=="T3"),], aes(y=perc.excess.CH4, x=plant.mass..g, 
                                                      fill=Treatment))+
  stat_smooth(aes(y=perc.excess.CH4,x=plant.mass..g, color=Treatment, linetype=Treatment), method="gam", alpha=0.2) +
  ylab(expression(paste("% excess CO"[2], sep=""))) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept= 100, linetype="longdash", color = "gray") +
  coord_cartesian(ylim=c(0, 1000)) +
  ggtitle("Day-59")+
  theme_classic()+
  theme(legend.position= "none")


#library(gridExtra)
line.excess<-plot_grid(T0.CO2, T1.CO2, T2.CO2, T3.CO2,
                       T0.CH4, T1.CH4, T2.CH4, T3.CH4, labels=c('A', '', '', '', 'B'), label_size=8, nrow=2, ncol=4)

line.excess
**Figure**. Percent excess concentrations for greenhouse gasses (**A**) carbon dioxide (CO~2~) and (**B**) methane (CH~4~) in tanks receiving burned and unburned plant material at the beginning of the experiment and during three experimental time points.  Gray line indicates 100% saturation under atmospheric equilibrium.

Figure. Percent excess concentrations for greenhouse gasses (A) carbon dioxide (CO2) and (B) methane (CH4) in tanks receiving burned and unburned plant material at the beginning of the experiment and during three experimental time points. Gray line indicates 100% saturation under atmospheric equilibrium.

ggsave("figures/GHG.excess.pdf", height=6, width=8)